Section 2 Scenario Analysis: WCS Priority Landscapes

# turn off the s2 processing 
## https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data
sf::sf_use_s2(FALSE)

See the USFS Wildfire Crisis Strategy (WCS) document for background information.

2.1 Scenarios Considered

There are three scenarios considered in this analysis of the constraints to mechanical forest health and fuel reduction treatments:

# mechanical mgmt allowed if:
n_sc_temp <- 3
data.frame(
  scenario = 1:3
  , cover = rep("Forest or Shrubland",n_sc_temp)
  , protected = rep("Not Protected",n_sc_temp)
  , slope = c("<40%","<60%","<60%")
  , road = c("<1,000ft","<2,000ft","<2,000ft")
  , riparian = c(">100ft",">100ft",">50ft")
  , admin = c("No Designation","No Designation","Any Designation")
) %>% 
  tidyr::pivot_longer(
    cols = -c(scenario)
  ) %>% 
  tidyr::pivot_wider(
    names_from = scenario
    , values_from = value
    , names_prefix = "scenario_"
  ) %>% 
  dplyr::mutate(
    name = factor(
      name
      , levels = c(
        "cover"
        , "protected"
        , "slope"
        , "road"
        , "riparian"
        , "admin"
      )
      , labels = c(
        "NLCD Cover Type"
        , "Protected or<br>IRA Status"
        , "Slope"
        , "Distance to<br>Nearest Road"
        , "Riparian<br>Buffer"
        , "Administrative<br>Designation"
      )
      , ordered = T
    )
  ) %>% 
  dplyr::arrange(name) %>% 
  dplyr::mutate(
    lvl = paste0("L",dplyr::row_number()-1,":")
  ) %>% 
  dplyr::relocate(lvl) %>% 
# make table
kableExtra::kable(
    caption = "Scenarios Considered for Determining Mechanical Treatment Operability<br>Mechanical treatment allowed if:"
    , col.names = c(
      ""
      , ""
      , "Scenario 1<br>Most Constrained"
      , "Scenario 2<br>Moderate"
      , "Scenario 3<br>Least Constrained"
    )
    , escape = F
  ) %>%  
  kable_classic(full_width=T) %>% 
  kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE)

Table 2.1: Scenarios Considered for Determining Mechanical Treatment Operability
Mechanical treatment allowed if:

		Scenario 1 Most Constrained	Scenario 2 Moderate	Scenario 3 Least Constrained
L0:	NLCD Cover Type	Forest or Shrubland	Forest or Shrubland	Forest or Shrubland
L1:	Protected or IRA Status	Not Protected	Not Protected	Not Protected
L2:	Slope	<40%	<60%	<60%
L3:	Distance to Nearest Road	<1,000ft	<2,000ft	<2,000ft
L4:	Riparian Buffer	>100ft	>100ft	>50ft
L5:	Administrative Designation	No Designation	No Designation	Any Designation

The figure below illustrates the workflow used to quantify the amount of land available for mechanical forest health and risk reduction fuel treatments by considering layered operational constraints. This example illustrates scenario 2 for a 10,252 hectare (25,000 acre) fireshed project area near Rustic, Colorado in the Colorado Front Range priority landscape.

2.2 Data Preparation

2.2.1 Load WCS landscapes

WCS priority landscape spatial data from the USFS Geospatial Data Discovery site (accessed 2023-04-04).

# load data
wf_landscapes <- sf::read_sf("../data/Wildfire_Crisis_Strategy_Landscapes/Wildfire_Crisis_Strategy_Landscapes_(Feature_Layer).shp") %>% 
  dplyr::rename_with(tolower) %>% 
  sf::st_make_valid() %>% 
  dplyr::mutate(
    hectares = (as.numeric(sf::st_area(geometry))/10000) 
    , Mil.Hectares = hectares %>% 
      scales::comma(suffix = " M", scale = 1e-6, accuracy = .01)
    , acres = (as.numeric(sf::st_area(geometry))/4046.85642) 
    , Mil.Acres = acres %>% 
      scales::comma(suffix = " M", scale = 1e-6, accuracy = .01)
  ) %>% 
  dplyr::left_join(
    data.frame(state = datasets::state.name, state_abb = datasets::state.abb)
    , by = join_by("state")
  ) %>% 
  dplyr::mutate(
    area_name = paste0(
      ifelse(is.na(state_abb) | state_abb == "", "UNK", state_abb)
      , ": "
      , name
    )
  ) %>% 
  # manual label placement
  dplyr::left_join(
    readr::read_csv("../data/wildfirepriority/area_label_placement.csv")
    , by = dplyr::join_by("area_name")
  )
#rename sf geom column
  names(wf_landscapes)[names(wf_landscapes)==tolower(attr(wf_landscapes, "sf_column"))] = "geometry"
  sf::st_geometry(wf_landscapes) = "geometry"
# set crs
  transform_crs <- sf::st_crs(wf_landscapes)
# view
  wf_landscapes %>% dplyr::glimpse()

2.2.2 Load landscape-level constraint data

Landscape-level constraint analysis data (i.e., row unique by WCS landscape) was created via this Google Earth Engine script.

# load landscape constraint scenario data
constrained_by_scnro_ls <- list.files("../data/wildfirepriority/all_slp/",pattern = "\\.csv$") %>%
  purrr::keep(stringr::str_starts(.,"wfpriority_all_slp_sc")) %>% 
  purrr::map(function(x){
    readr::read_csv(
      paste0("../data/wildfirepriority/all_slp/",x)
      , name_repair = "universal"
      , col_types = cols(.default = "c")
    ) %>% 
    dplyr::rename_with(tolower) %>% 
    dplyr::rename_with(make.names) %>% 
    dplyr::select(state,name,tidyselect::ends_with("_m2"),tidyselect::starts_with("pct_")) %>% 
    dplyr::mutate(scenario_id = stringr::word(x,4,sep="_") %>% readr::parse_number()) %>% 
    dplyr::relocate(scenario_id)
  }) %>%
  dplyr::bind_rows() %>%
  dplyr::inner_join(
    wf_landscapes %>%
      sf::st_drop_geometry() %>%
      dplyr::select(state,name,area_name)
    , by = dplyr::join_by(state,name)
  ) %>%
  dplyr::relocate(area_name,.after = "scenario_id") %>%
  dplyr::group_by(area_name,scenario_id) %>%
  dplyr::filter(dplyr::row_number()==1) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(
    dplyr::across(
      tidyselect::ends_with("_m2")
      , ~ as.numeric(.x) / 10000
    )
    , dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ as.numeric(.x)
    )
  ) %>%
  dplyr::rename_with(
    ~ gsub("_m2", "_ha", .x)
    , tidyselect::ends_with("_m2")
  ) %>%
  # calculate pct reduction
  dplyr::mutate(
    pct_rdctn1_protected = -1*(1 - pct_rmn1_protected)
    , pct_rdctn2_slope = -1*(pct_rmn1_protected - pct_rmn2_slope)
    , pct_rdctn3_roads = -1*(pct_rmn2_slope - pct_rmn3_roads)
    , pct_rdctn4_riparian = -1*(pct_rmn3_roads - pct_rmn4_riparian)
    , pct_rdctn5_administrative = -1*(pct_rmn4_riparian - pct_rmn5_administrative)
    , pct_rdctn_total = -1*(1 - pct_rmn5_administrative)
    , pct_covertype_area = covertype_area_ha/feature_area_ha
    # scenario description
    , scenario_lab = factor(
      scenario_id
      , levels = 1:3
      , labels = paste0("Scenario ", 1:3)
      , ordered = T
    ) %>% forcats::fct_rev()
    , scenario_desc = factor(
      scenario_id
      , levels = 1:3
      , labels = c("Scenario 1\n(status quo)", "Scenario 2\n(slopes+roads)", "Scenario 3\n(slopes+roads+admin)")
      , ordered = T
    )
  )

2.2.3 Load fireshed project area data

The fireshed registry spatial data was obtained from the USFS Geospatial Data Discovery tool for firesheds and project areas (accessed 2023-05-03).

See this section for exploration of fireshed and project area spatial data

### fireshed spatial data
fireshed <- sf::st_read("../data/firesheds/Fireshed_Registry3A_Fireshed/Fireshed_Registry%3A_Fireshed_(Feature_Layer).shp") %>%
  sf::st_transform(transform_crs) %>% 
  setNames(c(
      "shape_id"
      , "area_ha"
      , "fireshed_id"
      , "fireshed_name"
      , "fireshed_code"
      , "fireshed_state"
      , "nopas"
      , "objectid"
      , "fshed_id"
      , "exp_total"
      , "exp_usfs"
      , "exp_nonfs"
      , "exp_usfs_protected"
      , "exp_nonfs_protected"
      , "exp_usfs_managed"
      , "exp_nonfs_managed"
      , "exp_usfs_forest"
      , "exp_nonfs_forest"
      , "exp_usfs_nonforest"
      , "exp_nonfs_nonforest"
      , "exp_usfs_conifer"
      , "exp_nonfs_conifer"
      , "exp_usfs_managedforest"
      , "exp_nonfs_managedforest"
      , "exp_usfs_managedconifer"
      , "exp_nonfs_managedconifer"
      , "exp_nonfs_nonconifer_hihaz"
      , "dist_vs"
      , "crisis_strategy"
      , "key_preformance_indicator"
      , "national_usfs_rank"
      , "national_all_land_rank"
      , "regional_usfs_rank"
      , "regional_all_land_rank"
      , "start_date"
      , "end_date"
      , "geometry"
  )) %>% 
  dplyr::mutate(
    exposure_pct_rank = dplyr::percent_rank(exp_total)
    , exposure_pct_rank_grp = dplyr::case_when(
      exposure_pct_rank >= 1-0.01 ~ "Top 1%"
      , exposure_pct_rank >= 1-0.05 ~ "Top 5%"
      , exposure_pct_rank >= 1-0.10 ~ "Top 10%"
      , exposure_pct_rank >= 1-0.25 ~ "Top 25%"
      , TRUE ~ "Bottom 75%"
    ) %>% 
    factor(
      levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
      , ordered = T
    )
    # there is also a national_all_land_rank column
    , ntllandrank_pct_rank = dplyr::percent_rank(-national_all_land_rank)
    , ntllandrank_pct_rank_grp = dplyr::case_when(
        ntllandrank_pct_rank >= 1-0.01 ~ "Top 1%"
        , ntllandrank_pct_rank >= 1-0.05 ~ "Top 5%"
        , ntllandrank_pct_rank >= 1-0.10 ~ "Top 10%"
        , ntllandrank_pct_rank >= 1-0.25 ~ "Top 25%"
        , TRUE ~ "Bottom 75%"
      ) %>% 
      factor(
        levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
        , ordered = T
      )
    , crisis_strategy = ifelse(is.na(crisis_strategy),"Not High Risk",crisis_strategy) %>% 
      as.factor() %>% 
      forcats::fct_shift()
  )
  #rename sf geom column
    names(fireshed)[names(fireshed)==tolower(attr(fireshed, "sf_column"))] = "geometry"
    sf::st_geometry(fireshed) = "geometry"
    # calculate area
    fireshed <- fireshed %>% 
      dplyr::mutate(
        fireshed_area_ha = as.numeric(sf::st_area(geometry))/10000
        , fireshed_area_acres = (fireshed_area_ha*10000)/4046.85642
      )
## fireshed_proj_area spatial data
fireshed_proj_area <- sf::st_read("../data/firesheds/Fireshed_Registry3A_Project_Area/Fireshed_Registry%3A_Project_Area_(Feature_Layer).shp") %>%
  sf::st_transform(transform_crs) %>% 
  setNames(c(
      "shape_id"
      , "fireshed_id"
      , "pa_id"
      , "pa_area_ha"
      , "objectid"
      , "pa_id2"
      , "fshed_id"
      , "exp_total"
      , "exp_usfs"
      , "exp_nonfs"
      , "exp_usfs_protected"
      , "exp_nonfs_protected"
      , "exp_usfs_managed"
      , "exp_nonfs_managed"
      , "exp_usfs_forest"
      , "exp_nonfs_forest"
      , "exp_usfs_nonforest"
      , "exp_nonfs_nonforest"
      , "exp_usfs_conifer"
      , "exp_nonfs_conifer"
      , "exp_usfs_managedforest"
      , "exp_nonfs_managedforest"
      , "exp_usfs_managedconifer"
      , "exp_nonfs_managedconifer"
      , "exp_nonfs_nonconifer_hihaz"
      , "dist_vs"
      , "pctrecentlydisturbed"
      , "start_date"
      , "end_date"
      , "geometry"
  )) %>% 
  dplyr::mutate(
    exposure_pct_rank = dplyr::percent_rank(exp_total)
    , exposure_pct_rank_grp = dplyr::case_when(
      exposure_pct_rank >= 1-0.01 ~ "Top 1%"
      , exposure_pct_rank >= 1-0.05 ~ "Top 5%"
      , exposure_pct_rank >= 1-0.10 ~ "Top 10%"
      , exposure_pct_rank >= 1-0.25 ~ "Top 25%"
      , TRUE ~ "Bottom 75%"
    ) %>% 
    factor(
      levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
      , ordered = T
    )
  )
  #rename sf geom column
    names(fireshed_proj_area)[names(fireshed_proj_area)==tolower(attr(fireshed_proj_area, "sf_column"))] = "geometry"
    sf::st_geometry(fireshed_proj_area) = "geometry"
    # calculate area
    fireshed_proj_area <- fireshed_proj_area %>% 
      dplyr::mutate(
        pa_area_ha = as.numeric(sf::st_area(geometry))/10000
        , pa_area_acres = (pa_area_ha*10000)/4046.85642
      ) %>% 
      # JOIN WITH FIRESHED DATA
      dplyr::inner_join(
        fireshed %>%
          sf::st_drop_geometry() %>%
          dplyr::select(fireshed_id, crisis_strategy, exp_total
                        , exposure_pct_rank, exposure_pct_rank_grp
          ) %>% 
          dplyr::rename(exposure_total=exp_total) %>% 
          dplyr::rename_with(
            ~ paste0("fireshed_",.x)
            , -c(fireshed_id)
          )
        , by = dplyr::join_by(fireshed_id)
      ) %>%
      dplyr::select(pa_id,pa_area_ha
                    ,exp_total,exposure_pct_rank,exposure_pct_rank_grp
                    , tidyselect::starts_with("fireshed_")
      ) %>% 
      dplyr::rename(exposure_total=exp_total) %>%
      dplyr::rename_with(
        ~ paste0("pa_", .x, recycle0 = TRUE)
        , tidyselect::starts_with("exp")
      )

Fireshed project area constraint analysis data (i.e., row unique by WCS landscape + f.p.a.) was created via this Google Earth Engine script

# load fireshed constraint scenario data
constrained_by_scnro_ls_pa <-
  list.files("../data/wildfirepriority/fireshed_slp/",pattern = "\\.csv$") %>%
  purrr::keep(stringr::str_starts(.,"wfpriority_fireshed_slp_sc")) %>% 
  purrr::map(function(x){
    readr::read_csv(
      paste0("../data/wildfirepriority/fireshed_slp/",x)
      , name_repair = "universal"
      , col_types = cols(.default = "c")
    ) %>% 
    dplyr::rename_with(tolower) %>% 
    dplyr::rename_with(make.names) %>% 
    dplyr::select(state,name,pa_id,tidyselect::ends_with("_m2"),tidyselect::starts_with("pct_")) %>% 
    dplyr::mutate(scenario_id = stringr::word(x,4,sep="_") %>% readr::parse_number()) %>% 
    dplyr::relocate(scenario_id)
  }) %>%
  dplyr::bind_rows() %>%
  dplyr::inner_join(
    wf_landscapes %>%
      sf::st_drop_geometry() %>%
      dplyr::select(state,name,area_name)
    , by = dplyr::join_by(state,name)
  ) %>%
  dplyr::relocate(area_name,.after = "scenario_id") %>%
  dplyr::mutate(
    dplyr::across(
      tidyselect::ends_with("_m2")
      , ~ as.numeric(.x) / 10000
    )
    , dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ as.numeric(.x)
    )
  ) %>%
  dplyr::rename_with(
    ~ gsub("_m2", "_ha", .x)
    , tidyselect::ends_with("_m2")
  ) %>%
  dplyr::group_by(area_name,pa_id,scenario_id) %>%
  dplyr::arrange(desc(pct_pa_intrsct)) %>%
  dplyr::filter(dplyr::row_number()==1) %>%
  dplyr::ungroup() %>%
  dplyr::filter(pct_pa_intrsct>=0.00001) %>%
  # calculate pct reduction
  dplyr::mutate(
    pct_rdctn1_protected = -1*(1 - pct_rmn1_protected)
    , pct_rdctn2_slope = -1*(pct_rmn1_protected - pct_rmn2_slope)
    , pct_rdctn3_roads = -1*(pct_rmn2_slope - pct_rmn3_roads)
    , pct_rdctn4_riparian = -1*(pct_rmn3_roads - pct_rmn4_riparian)
    , pct_rdctn5_administrative = -1*(pct_rmn4_riparian - pct_rmn5_administrative)
    , pct_rdctn_total = -1*(1 - pct_rmn5_administrative)
    , pct_covertype_area = covertype_area_ha/feature_area_ha
    # calculate level of constraint
      , cnstrnt_lvl = dplyr::case_when(
          -pct_rdctn_total >= 0.85 ~ 1
          , -pct_rdctn_total >= 0.65 ~ 2
          , -pct_rdctn_total >= 0.0 ~ 3
        )
      , cnstrnt_class = factor(
          cnstrnt_lvl 
          , levels = 1:3
          , labels = c("high constraint", "med. constraint", "low constraint")
          , ordered = T
        ) %>% forcats::fct_rev()
      , rmn_cnstrnt_class = factor(
          cnstrnt_lvl 
          , levels = 1:3
          , labels = c("0–15% treatable", "16–35% treatable", ">35% treatable")
          , ordered = T
        ) %>% forcats::fct_rev()
      # scenario description
      , scenario_lab = factor(
          scenario_id
          , levels = 1:3
          , labels = paste0("Scenario ", 1:3)
          , ordered = T
        ) %>% forcats::fct_rev()
      , scenario_desc = factor(
        scenario_id
        , levels = 1:3
        , labels = c("Scenario 1\n(status quo)", "Scenario 2\n(slopes+roads)", "Scenario 3\n(slopes+roads+admin)")
        , ordered = T
      )
  ) %>% 
  # join with fireshed data
  dplyr::inner_join(
    fireshed_proj_area %>%
      sf::st_drop_geometry() %>%
      dplyr::select(pa_id, tidyselect::starts_with("pa_exposure"), tidyselect::starts_with("fireshed_")) %>%
      dplyr::mutate(pa_id=as.character(pa_id))
    , by = dplyr::join_by(pa_id)
  )

2.3 Geographic Delineation Description

Timeline of the research utilized to motivate and support the creation of the WCS:

• 2018 August - Toward share stewardship across the landscapes: an outcome-based investment strategy created a hotspot map of the source of wildfire risk to communities showing that 80% of community exposure from NFS lands originates on 20% of the NFS area.
• 2019 February - Evers et al. analyzed community wildfire exposure from national forests by associating conditions that affect exposure in the areas where wildfires ignite to conditions where exposure likely occurs
• 2019 May - Ager et al. mapped community firesheds by creating a continuous smoothed surface of predicted structure exposure from all FSim ignitions (see also Palaiologou et al. 2019 and Ager et al. 2019)
• 2020 - Evers et al. release Fireshed Registry as a geospatial dashboard and decision tool
• 2021 May - Ager et al. 2021 official release of the Fireshed Registry as a geospatial information system that organizes landscape risk to developed areas into a planning framework (see figure below)
• 2022 January - US Department of Agriculture, Forest Service launches a 10-year Wildfire Crisis Strategy
• 2022 April - 10 initial landscapes selected for WCS prioritization
• 2023 January - additional 11 landscapes selected for WCS prioritization

Nested spatial framework for firesheds. Each scale has specific functionality in terms of the planning processes. Firesheds are the broad scale unit of prioritization, but project areas within them are also prioritized as part of implementation of treatments. Project areas are roughly the size that national forests use for conducting vegetation and fuel management projects. The relative variation among firesheds compared to variation within them controls the relative emphasis on selecting firesheds versus individual planning areas. Ager et al. 2021

#define basemap
states_temp <- USAboundaries::us_states(
  states = c(
        # usfs region 1-6 states
        "MT","WY","CO","NM","AZ","UT","ID","WA","OR","CA","NV"
        # , "KS","NE","SD","ND"
      )
  ) %>% 
  sf::st_transform(transform_crs)
cities_temp <- USAboundaries::us_cities(
  states = c(
        # usfs region 1-6 states
        "MT","WY","CO","NM","AZ","UT","ID","WA","OR","CA","NV"
        # , "KS","NE","SD","ND"
      )
  ) %>% 
  sf::st_transform(transform_crs) %>% 
  dplyr::filter(
    toupper(city) %in% c(
      "DENVER"
      , "TUSCON"
      , "SALT LAKE CITY"
      , "LAS VEGAS"
      , "SAN DIEGO"
      , "LOS ANGELES"
      , "FRESNO"
      , "SAN FRANCISCO"
      , "SACRAMENTO"
      , "PORTLAND"
      , "SEATTLE"
      , "ALBUQUERQUE"
      , "RENO"
      , "BOISE"
      , "HELENA"
      , "BILLINGS"
    )
    | (toupper(city) == "PHOENIX" & state_abbr == "AZ")
  )
# plot basemap
plt_basemap <-
ggplot() + 
  geom_sf(
    data = states_temp
    , fill = NA
    , color = "gray66"
  ) +
  geom_sf_text(
    data = states_temp
    , mapping = aes(label = stusps)
    , color = "gray66"
    , size = 3
  ) +
  geom_sf(
    data = cities_temp
    , shape = 1
    , size = 1
    , color = "gray44"
  ) +
  geom_sf_text(
    data = cities_temp
    , mapping = aes(label = city)
    , color = "gray44"
    , size = 2
    # , hjust = -0.15
    , hjust = 0.5
    , vjust = -0.25
  ) +
  coord_sf(expand = F) +
  theme_void()

# plot basemap simple
plt_basemap_simple <-
ggplot() + 
  geom_sf(
    data = states_temp
    , fill = NA
    , color = "gray66"
  ) +
  coord_sf(expand = F) +
  theme_void()

2.3.1 Fireshed risk to communities (2021)

ggplot() + 
  geom_sf(
    data = fireshed
    , mapping = aes(fill=ntllandrank_pct_rank_grp)
  ) + 
  geom_sf(
    data = USAboundaries::us_states() %>% 
      dplyr::filter(!stusps %in% c("AK","HI","PR")) %>% 
      sf::st_transform(transform_crs)
    , fill = NA
    , color = "gray88"
  ) +
  scale_fill_manual(values=c("firebrick","orange3","khaki","seagreen3","royalblue")) +
  coord_sf(expand = F) +
  labs(
    fill = "Fireshed exposure ranks"
  ) +
  theme_void() + 
  theme(
    legend.position = "top"
    , legend.direction = "horizontal"
    , legend.title = element_text(size = 8)
    , legend.text = element_text(size = 7)
  )

National map of the 7,688 firesheds created from community wildfire transmission data (Evers et al. 2020). The fireshed boundaries were created with a process that delineates hotspots of fire transmission to buildings in adjacent or nearby communities. See the Methods section for details on delineating firesheds. Figure 2 from Ager et al. 2021 (p. 7)

2.3.2 Wilfire Crisis Strategy 21 Priority Landscapes (2022-2023)

plt_fshed <- 
plt_basemap + 
  geom_sf(
    data = fireshed %>%
      sf::st_filter(
        wf_landscapes %>%
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(fill=crisis_strategy)
    , alpha = 0.8
  ) +
  geom_sf(
    data = wf_landscapes
    , mapping = aes(color = paste0(investment, " investment"))
    , fill=NA
    # , color="black" 
    , lwd=0.6
  )+
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = stringr::str_wrap(
          dplyr::case_when(
            stringr::str_count(name, "\\w+")<2
              ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
            , stringr::str_starts(name,"Sierra and Elko") ~ stringr::word(name,1,3)
            , TRUE ~ stringr::word(name,1,2)
          )
          , 7)
        , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , stat = "sf_coordinates"
    , size = 2.5
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
  scale_color_manual(values = c("black", "gray33")) +
  labs(
    fill = "Fireshed\nHigh-Risk Status"
    , color = "Landscapes"
  ) +
  theme(
      legend.position = c(0.87,0.87) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 8)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(lwd = 3, fill = NA)) #, linetype = c(5,1,1)
    , fill = guide_legend(order = 2, override.aes = list(lwd = 0))
  )
# plot
plt_fshed

The boundaries of the 21 priority landscapes roughly follow the boundaries of “firesheds” prioritized to reduce wildfire transmission to developed areas (Figure 1). Firesheds are geographic delineations with an average area of approximately 250,000 acres (~100,000 hectares) that were created to organize the landscape into units for managing wildfire risk to communities.

2.3.3 Wilfire Crisis Strategy Treatment Objective

plt_basemap + 
  geom_sf(
    data = fireshed %>%
      dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
      sf::st_filter(
        wf_landscapes %>%
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(fill=crisis_strategy)
    , alpha = 0.8
  ) +
  geom_sf(
    data = wf_landscapes
    , mapping = aes(color = paste0(investment, " investment"))
    , fill=NA
    # , color="black" 
    , lwd=0.6
  )+
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = stringr::str_wrap(
          dplyr::case_when(
            stringr::str_count(name, "\\w+")<2
              ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
            , stringr::str_starts(name,"Sierra and Elko") ~ stringr::word(name,1,3)
            , TRUE ~ stringr::word(name,1,2)
          )
          , 7)
        , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , stat = "sf_coordinates"
    , size = 2.5
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
  scale_color_manual(values = c("black", "gray33")) +
  labs(
    fill = "Fireshed\nHigh-Risk Status"
    , color = "Landscapes"
  ) +
  theme(
      legend.position = c(0.87,0.87) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 8)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(lwd = 3, fill = NA)) #, linetype = c(5,1,1)
    , fill = guide_legend(order = 2, override.aes = list(lwd = 0))
  )

Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the Strategy (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).

2.3.4 Fireshed Project Areas

name_temp <- "Enchanted Circle"
# pa map
plt_fshed_pa <-
ggplot() + 
  geom_sf(data = 
    fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(
            # fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
            name == name_temp
          ) %>% 
          dplyr::select(pa_id)
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      )
    , aes(fill = fireshed_crisis_strategy, color = "Project Area\nboundary")
    , alpha = 0.8
    , lwd = 1.2
  ) +
  geom_sf(
    data = fireshed %>%
      # dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
      sf::st_filter(
        wf_landscapes %>%
          dplyr::filter(name == name_temp) %>% 
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(color = "Fireshed\nboundary")
    , fill = NA
    , lwd = 1.2
    , linetype = "dashed"
  ) +
  geom_sf(
    data = wf_landscapes %>% dplyr::filter(name == name_temp)
    , mapping = aes(color = "WCS Landscape\nboundary")
    , fill=NA
    , lwd = 0.9
  )+
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
  scale_color_manual(values = c("gray44","skyblue2", "black")) +
  labs(
    fill = "Fireshed\nHigh-Risk Status"
    , color = "Spatial\nFramework"
    , subtitle = name_temp
  ) +
  coord_sf(expand = F) +
  theme_void() +
  theme(
      legend.position = c(1.08,0.79) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 10)
      , plot.subtitle = element_text(face = "bold.italic", size = 10, hjust = 0.5)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(size = 8, linewidth = 5, fill = NA)) #, linetype = c(5,1,1)))
    , fill = guide_legend(order = 2, override.aes = list(size = 8,linewidth = 0, alpha = 1)) # , color = NA
  )

# pa map w/in
plt_fshed_pa_incl <- ggplot() + 
  geom_sf(data = 
    fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(
            # fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
            name == name_temp
          ) %>%
          dplyr::mutate(
            # is_25in = ifelse(pct_pa_intrsct>=0.25,"≥25%","<25%")
            is_25in = ifelse(pct_pa_intrsct>=0.25,"included","excluded")
          ) %>% 
          dplyr::select(
            pa_id
            , is_25in
          )
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      )
    , aes(fill = is_25in, color = "Project Area\nboundary")
    , alpha = 0.9
    , lwd = 1.2
  ) +
  geom_sf(
    data = fireshed %>%
      # dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
      sf::st_filter(
        wf_landscapes %>%
          dplyr::filter(name == name_temp) %>% 
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(color = "Fireshed\nboundary")
    , fill = NA
    , lwd = 1.2
    , linetype = "dashed"
  ) +
  geom_sf(
    data = wf_landscapes %>% dplyr::filter(name == name_temp)
    , mapping = aes(color = "WCS Landscape\nboundary")
    , fill=NA
    , lwd = 0.9
  )+
  scale_fill_manual(values = viridis::viridis(n=4, direction = -1, alpha = 0.9)[c(1,3)]) +
  scale_color_manual(values = c("gray44","skyblue2", "black")) +
  labs(
    fill = "Project Area\nInclusions"
    # fill = "Included in this analysis\nif ≥25% area w/in landscape"
    , color = "Spatial\nFramework"
    , subtitle = name_temp
  ) +
  coord_sf(expand = F) +
  theme_void() +
  theme(
      legend.position = c(1.08,0.83) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 10)
      , plot.subtitle = element_text(face = "bold.italic", size = 10, hjust = 0.5)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(size = 8, linewidth = 5, fill = NA)) #, linetype = c(5,1,1)))
    , fill = guide_legend(order = 2, override.aes = list(size = 8,linewidth = 0, alpha = 1)) # , color = NA
  )

# combine
cowplot::plot_grid(
    plotlist = list(
      NULL
      , plt_fshed_pa + 
        theme(
          plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
          , plot.margin = margin(0.5,3,0.5,0.5, "cm")
        )
      , NULL
      , plt_fshed_pa_incl + 
        theme(
          plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
          , plot.margin = margin(0.5,3,0.5,0.5, "cm")
        )
      , NULL
    )
    , nrow = 1
    , rel_widths = c(0.01,1,0.01,1,0.01)
    , labels = c("","A","","B","")
    , label_size = 17
    , label_fontface = "bold"
    , label_colour = "gray11"
    # , label_y = 0.99
  ) +
  theme(plot.background = element_rect(fill = "white", color = NA))

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_ex_spatial_strctr.png")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 5.8
    , units = "in"
    , dpi = "print"
  )

1. Nested within firesheds are project areas of approximately 10,000 ha in size which represent the geographic unit at which vegetation and fuel management projects are planned (Ager et al. 2021).

2. Fireshed project areas that did not have at least 25% of area within the boundary of the Strategy priority landscape were excluded from this analysis. This cutoff was implemented under the assumption that with less than 25% of area available, treatment alone would not substantially affect wildfire behavior across the fireshed (North et al. 2015). For each fireshed project area included in this analysis, the entire area within the project area boundary was used to calculate the percent of mechanically available acreage including area outside of the priority landscape boundary.

2.4 Landscape-level analysis

# filter fireshed project area
constrained_by_scnro_ls_pa <- constrained_by_scnro_ls_pa %>% 
  dplyr::filter(pct_pa_intrsct>=0.25)

2.4.1 Reduction Treatable Area Table

# scnum_temp<-2
  tbl_temp <- constrained_by_scnro_ls %>% 
    # dplyr::filter(scenario_id==scnum_temp) %>% 
    dplyr::mutate(
      dplyr::across(
        tidyselect::ends_with("_ha")
        , ~ scales::comma(.x, accuracy = 1)
      )
      , dplyr::across(
        tidyselect::starts_with("pct_")
        , ~ scales::percent(.x, accuracy = 0.1)
      )
    ) %>% 
    dplyr::select(
      area_name
      , scenario_lab
      , covertype_area_ha
      , pct_covertype_area
      , pct_rdctn1_protected
      , pct_rdctn2_slope
      , pct_rdctn3_roads
      , pct_rdctn4_riparian
      , pct_rdctn5_administrative
      , rmn5_administrative_area_ha
      , pct_rmn5_administrative
    ) %>% 
    dplyr::arrange(area_name,desc(scenario_lab))
  # make table 
  kableExtra::kable(
      tbl_temp %>% dplyr::select(-c(area_name))
      # , caption = paste0("<b><font color=navy>Scenario "
      #                    ,scnum_temp
      #                    ,"</font></b><br>percent reduction of different types of constraints on mechanical treatment"
      #         )
      , caption = "<b>Wildfire Crisis Strategy 21 Priority Landscapes</b><br>percent reduction of different types of constraints on mechanical treatment"
      , col.names = c(
        ""
        , "Forest+Shrub\n(ha)", "Forest+Shrub\n%"
        , "Protected"
        , "Slope\nSteepness"
        , "Road\nDistance"
        , "Riparian\nBuffer"
        , "Administrative"
        , "Remaining (ha)"
        , "Remaining (%)"
      )
      , escape = F
    ) %>%  
    add_header_above(c(" " = 1, "Area of\nPriority Landscape"=2, "Constraint\nLeast Flexible to Most Flexible" = 5, " " = 2)) %>%
    kable_classic(full_width=T) %>% 
    pack_rows(index = table(forcats::fct_inorder(tbl_temp$area_name))) %>% 
    kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%  
    kableExtra::scroll_box(width = "740px")

Table 2.2: Wildfire Crisis Strategy 21 Priority Landscapes
percent reduction of different types of constraints on mechanical treatment

	Area of Priority Landscape		Constraint Least Flexible to Most Flexible
	Forest+Shrub (ha)	Forest+Shrub %	Protected	Slope Steepness	Road Distance	Riparian Buffer	Administrative	Remaining (ha)	Remaining (%)
AZ: 4FRI
Scenario 1	2,190,804	90.6%	-13.8%	-7.5%	-26.4%	-10.0%	-15.2%	592,541	27.0%
Scenario 2	2,190,804	90.6%	-13.8%	-2.0%	-13.2%	-13.9%	-19.7%	819,126	37.4%
Scenario 3	2,190,804	90.6%	-13.8%	-2.0%	-13.2%	-9.8%	0.0%	1,339,032	61.1%
AZ: Prescott
Scenario 1	219,110	86.4%	-12.1%	-16.6%	-24.0%	-12.1%	-5.3%	65,793	30.0%
Scenario 2	219,110	86.4%	-12.1%	-3.2%	-12.2%	-18.9%	-8.3%	99,404	45.4%
Scenario 3	219,110	86.4%	-12.1%	-3.2%	-12.2%	-13.4%	0.0%	129,550	59.1%
AZ: San Carlos Apache Tribal Forest Protection
Scenario 1	1,056,065	91.4%	-6.1%	-20.4%	-51.9%	-5.8%	-0.8%	157,750	14.9%
Scenario 2	1,056,065	91.4%	-6.1%	-5.3%	-48.8%	-10.2%	-1.7%	294,497	27.9%
Scenario 3	1,056,065	91.4%	-6.1%	-5.3%	-48.8%	-7.3%	0.0%	342,834	32.5%
CA: North Yuba
Scenario 1	139,202	96.5%	-11.6%	-33.4%	-7.8%	-14.7%	-2.3%	42,013	30.2%
Scenario 2	139,202	96.5%	-11.6%	-13.9%	-3.2%	-23.0%	-3.3%	62,696	45.0%
Scenario 3	139,202	96.5%	-11.6%	-13.9%	-3.2%	-16.7%	0.0%	76,029	54.6%
CA: Plumas Community Protection
Scenario 1	100,872	93.3%	0.0%	-23.2%	-19.8%	-21.5%	-1.9%	33,888	33.6%
Scenario 2	100,872	93.3%	0.0%	-5.2%	-8.1%	-33.9%	-2.9%	50,328	49.9%
Scenario 3	100,872	93.3%	0.0%	-5.2%	-8.1%	-24.4%	0.0%	62,810	62.3%
CA: Southern California Fireshed Risk Reduction Strategy
Scenario 1	1,408,168	86.0%	-57.0%	-21.0%	-8.3%	-2.3%	-1.4%	139,165	9.9%
Scenario 2	1,408,168	86.0%	-57.0%	-9.6%	-6.5%	-4.2%	-3.0%	277,540	19.7%
Scenario 3	1,408,168	86.0%	-57.0%	-9.6%	-6.5%	-3.0%	0.0%	336,002	23.9%
CA: Stanislaus
Scenario 1	116,139	94.3%	-2.9%	-26.7%	-10.6%	-24.7%	-0.6%	40,057	34.5%
Scenario 2	116,139	94.3%	-2.9%	-10.7%	-3.9%	-35.8%	-0.8%	53,380	46.0%
Scenario 3	116,139	94.3%	-2.9%	-10.7%	-3.9%	-26.1%	0.0%	65,599	56.5%
CA: Trinity Forest Health and Fire-Resilient Rural Communities
Scenario 1	295,916	80.4%	-12.9%	-39.2%	-10.8%	-10.7%	-12.0%	42,598	14.4%
Scenario 2	295,916	80.4%	-12.9%	-11.9%	-6.5%	-20.7%	-24.1%	70,790	23.9%
Scenario 3	295,916	80.4%	-12.9%	-11.9%	-6.5%	-15.1%	0.0%	158,737	53.6%
CO: Colorado Front Range
Scenario 1	1,051,540	72.6%	-20.4%	-19.9%	-20.0%	-11.3%	-1.6%	281,048	26.7%
Scenario 2	1,051,540	72.6%	-20.4%	-5.6%	-12.0%	-17.3%	-3.2%	435,578	41.4%
Scenario 3	1,051,540	72.6%	-20.4%	-5.6%	-12.0%	-12.3%	0.0%	521,905	49.6%
ID: Nez Perce-Clearwater-Lower Salmon
Scenario 1	699,884	89.0%	-43.0%	-19.7%	-9.1%	-1.4%	-0.3%	185,161	26.5%
Scenario 2	699,884	89.0%	-43.0%	-6.7%	-5.5%	-2.3%	-0.7%	292,841	41.8%
Scenario 3	699,884	89.0%	-43.0%	-6.7%	-5.5%	-1.6%	0.0%	302,267	43.2%
ID: Southwest Idaho
Scenario 1	608,338	87.6%	-13.3%	-28.3%	-10.4%	-2.0%	-0.3%	278,231	45.7%
Scenario 2	608,338	87.6%	-13.3%	-7.0%	-6.6%	-2.8%	-0.8%	422,974	69.5%
Scenario 3	608,338	87.6%	-13.3%	-7.0%	-6.6%	-2.0%	0.0%	432,970	71.2%
MT: Kootenai Complex
Scenario 1	372,345	91.4%	-21.4%	-15.0%	-7.8%	-6.9%	-10.6%	142,439	38.3%
Scenario 2	372,345	91.4%	-21.4%	-3.6%	-2.1%	-8.6%	-13.8%	188,288	50.6%
Scenario 3	372,345	91.4%	-21.4%	-3.6%	-2.1%	-6.1%	0.0%	249,151	66.9%
NM: Enchanted Circle
Scenario 1	506,460	85.6%	-8.5%	-21.5%	-33.4%	-7.5%	-2.2%	135,772	26.8%
Scenario 2	506,460	85.6%	-8.5%	-5.2%	-25.7%	-11.7%	-3.7%	228,840	45.2%
Scenario 3	506,460	85.6%	-8.5%	-5.2%	-25.7%	-8.2%	0.0%	265,124	52.3%
NV: Sierra and Elko Fronts
Scenario 1	1,003,525	73.2%	-26.6%	-9.8%	-26.9%	-6.2%	-0.9%	297,023	29.6%
Scenario 2	1,003,525	73.2%	-26.6%	-2.0%	-15.8%	-8.9%	-1.5%	454,449	45.3%
Scenario 3	1,003,525	73.2%	-26.6%	-2.0%	-15.8%	-6.3%	0.0%	495,822	49.4%
OR: Central Oregon
Scenario 1	937,977	88.8%	-7.2%	-2.7%	-14.9%	-5.1%	-10.2%	562,234	59.9%
Scenario 2	937,977	88.8%	-7.2%	-0.5%	-4.7%	-6.0%	-11.6%	656,297	70.0%
Scenario 3	937,977	88.8%	-7.2%	-0.5%	-4.7%	-4.2%	0.0%	781,760	83.3%
OR: Klamath River Basin
Scenario 1	2,900,530	77.7%	-20.5%	-16.9%	-17.0%	-5.5%	-4.6%	1,032,327	35.6%
Scenario 2	2,900,530	77.7%	-20.5%	-5.7%	-8.1%	-8.6%	-8.2%	1,420,954	49.0%
Scenario 3	2,900,530	77.7%	-20.5%	-5.7%	-8.1%	-6.2%	0.0%	1,729,729	59.6%
OR: Mount Hood Forest Health and Fire-Resilient Communities
Scenario 1	334,397	76.7%	-19.6%	-14.2%	-16.9%	-9.1%	-16.4%	79,550	23.8%
Scenario 2	334,397	76.7%	-19.6%	-4.1%	-7.3%	-13.2%	-20.9%	116,660	34.9%
Scenario 3	334,397	76.7%	-19.6%	-4.1%	-7.3%	-9.5%	0.0%	198,932	59.5%
UT: Pine Valley
Scenario 1	153,246	94.2%	-37.2%	-6.9%	-21.5%	-6.7%	0.0%	42,532	27.8%
Scenario 2	153,246	94.2%	-37.2%	-1.4%	-10.5%	-9.2%	0.0%	63,919	41.7%
Scenario 3	153,246	94.2%	-37.2%	-1.4%	-10.5%	-6.5%	0.0%	68,040	44.4%
UT: Wasatch
Scenario 1	381,429	89.5%	-43.0%	-12.1%	-14.8%	-4.9%	-0.3%	95,232	25.0%
Scenario 2	381,429	89.5%	-43.0%	-3.7%	-7.9%	-6.7%	-0.5%	145,982	38.3%
Scenario 3	381,429	89.5%	-43.0%	-3.7%	-7.9%	-4.8%	0.0%	155,181	40.7%
WA: Central Washington Initiative
Scenario 1	670,440	68.1%	-20.6%	-32.0%	-12.3%	-8.5%	-12.2%	96,766	14.4%
Scenario 2	670,440	68.1%	-20.6%	-10.4%	-7.8%	-14.1%	-21.8%	168,986	25.2%
Scenario 3	670,440	68.1%	-20.6%	-10.4%	-7.8%	-10.2%	0.0%	341,540	50.9%
WA: Colville Northeast Washington Vision
Scenario 1	655,989	92.8%	-20.1%	-18.5%	-16.0%	-10.1%	-1.2%	223,368	34.1%
Scenario 2	655,989	92.8%	-20.1%	-3.8%	-7.5%	-14.4%	-2.0%	342,407	52.2%
Scenario 3	655,989	92.8%	-20.1%	-3.8%	-7.5%	-10.4%	0.0%	382,335	58.3%

2.4.2 Change in Available by Scenario

constrained_by_scnro_ls %>% 
  dplyr::mutate(scenario_lab = scenario_lab %>% forcats::fct_rev()) %>% 
ggplot(
    mapping = aes(x=scenario_lab,y=pct_rmn5_administrative,group=area_name)
    ) +
    geom_line(mapping=aes(color = area_name), linewidth = 1.5) +
    geom_label(
      mapping=aes(label = scales::percent(pct_rmn5_administrative, scale = 100, accuracy = 1))
      , color = "black"
      , size = 3
      , label.padding = unit(0.15, "lines")
    ) +
    facet_wrap(facets = vars(area_name)
      , ncol = 3
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
      , scales = "free_x"
    ) +
    scale_color_viridis_d(option = "turbo", alpha = 0.8) +
    scale_y_continuous(expand = expansion(mult = c(0.1, 0.1)), labels = scales::percent_format(), breaks = scales::extended_breaks(6)) +
    # scale_x_discrete(labels = scales::label_wrap(20)) +
    labs(
      x = "Most Restrictive \U2192 Least Restrictive"
      , y = "Percent Treatable Area Remaining"
    ) +
    theme_light() +
    theme(
      legend.position = "none"
      , axis.text.y = element_text(size=7)
      , axis.text.x = element_text(size=7.5)
      , strip.text = element_text(color = "black", face = "bold", size = 10)
    )

2.4.3 Reduction by Constraint

# reshape
constrained_by_scnro_ls_long <- constrained_by_scnro_ls %>%
  dplyr::select(state, name, area_name, tidyselect::starts_with("scenario_"), tidyselect::starts_with("pct_rdctn")) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("pct_rdctn")
    , names_to = "constraint"
    , values_to = "pct_rdctn"
    , names_prefix = "pct_rdctn"
    , values_drop_na = F
  ) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , constraint = factor(
        constraint
        , ordered = TRUE
        , levels = c(
          "protected"
          , "slope"
          , "roads"
          , "riparian"
          , "administrative"
          , "total"   
        )
        , labels = c(
          "Protected"
          , "Slope\nSteepness"
          , "Road\nDistance"
          , "Riparian\nBuffer"
          , "Administrative"
          , "Total"
        )
      ) %>% forcats::fct_rev()
  ) %>% 
  dplyr::left_join(
    constrained_by_scnro_ls %>% dplyr::select(state,name,scenario_id,pct_rdctn_total)
    , by = join_by(state,name,scenario_id)
  )

# color pallete removing dark blue
n_temp <- ((constrained_by_scnro_ls_long$constraint %>% unique() %>% length())-1)*2
plasma_custom <- viridis::plasma(n = n_temp)[seq(2,n_temp,2)]
# scales::show_col(plasma_custom)
# plasma_custom
# plot
ggplot(
  data = constrained_by_scnro_ls_long %>% 
      dplyr::filter(constraint!="Total") %>% 
      dplyr::mutate(
        dplyr::across(
          c(scenario_lab, area_name)
          , ~ forcats::fct_rev(.x)
        )
      )
  , mapping = aes(y = area_name, x = pct_rdctn)
) +
  geom_col(
    mapping = aes(fill = constraint)
    , color = NA, width = 0.7
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      label = scales::percent(ifelse(pct_rdctn<0.11*-1,pct_rdctn,NA), accuracy = 1)
      , group = constraint
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_long %>% 
      dplyr::filter(constraint=="Total") %>% 
      dplyr::mutate(
        dplyr::across(
          c(scenario_lab, area_name)
          , ~ forcats::fct_rev(.x)
        )
      )
    , mapping = aes(
      y = area_name, x = pct_rdctn_total
      , label = scales::percent(pct_rdctn_total, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_grid(cols = vars(scenario_lab)) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0.04, .3)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x =  element_blank() #
    , strip.text = element_text(color = "black", face = "bold")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

ggplot2::ggsave(
    filename = paste0("../data/plt_cnstrnt_rdctn_ls.png")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

2.5 Fireshed-level analysis

Based on model simulations of how much area generally needs to be treated to influence wildfire behavior, we binned the firesheds into three classes of mechanical constraint:

• high (85–100% [i.e., only 0–15% is available for mechanical treatment]): fuels treatment would principally need to rely on fire
• medium (65– 84%): could use a combination of fire and mechanical thinning
• low (65%): could effectively influence wildfire behavior with mechanical treatment alone

2.5.1 Distribution of Fireshed PA Constraint

constrained_by_scnro_ls_pa %>% 
  dplyr::count(area_name,scenario_lab,cnstrnt_class) %>% 
  dplyr::group_by(area_name,scenario_lab) %>% 
  dplyr::mutate(
    pct = n/sum(n)
    , high_pct=max(ifelse(cnstrnt_class=="high constraint",pct,0)) 
    , tot = sum(n)
    , scenario_lab = scenario_lab %>% forcats::fct_rev()
  ) %>% 
ggplot() +
  geom_col(
    mapping = aes(x = pct, y = reorder(area_name, desc(area_name)), fill=cnstrnt_class)
    ,width = 0.7, alpha=0.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = reorder(area_name, desc(area_name)), group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.055,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  # geom_text(
  #     data = constrained_byftr_huc12_wide %>%
  #       dplyr::group_by(area_name) %>%
  #       dplyr::summarise(n=n(), high_n=sum(ifelse(cnstrnt_class=="high constraint",1,0)))
  #     , mapping = aes(x = -0.05,y=reorder(area_name, -(high_n/n)),label=paste0("n=",scales::comma(n)))
  #     , size = 3
  #     , color = "black"
  #     , hjust = 0.7
  #   ) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  facet_grid(cols = vars(scenario_lab)) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "% of Fireshed\nProject Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

2.5.2 Distribution of Fireshed Area

# plot
ggplot(data = constrained_by_scnro_ls_pa %>% dplyr::filter(scenario_id==1)
         , mapping = aes(x = feature_area_ha, y = reorder(area_name,desc(area_name) ))
    ) +
    geom_vline(
      xintercept = 
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(scenario_id==1) %>% 
          dplyr::pull(feature_area_ha) %>% 
          median()
      , linetype="dashed", color="gray66"
    ) +
    geom_violin(aes(fill = state)) + 
    geom_boxplot(width = 0.15) +
    # geom_point(size = 0.7, color = "gray33", alpha = 0.2) +
    geom_text(
      data = constrained_by_scnro_ls_pa %>% 
        dplyr::filter(scenario_id==1) %>% 
        dplyr::group_by(area_name) %>% 
        dplyr::summarise(n=n(), max_area=max(feature_area_ha))
      , mapping = aes(
          x = (
            constrained_by_scnro_ls_pa %>% 
            dplyr::filter(scenario_id==1) %>% 
            dplyr::pull(feature_area_ha) %>% 
            min()
          )*.9
          ,y=reorder(area_name,desc(area_name) )
          ,label=paste0("n=",scales::comma(n))
        )
      , size = 3
      , color = "black"
      , hjust = 0.7
    ) +
    scale_fill_viridis_d(option = "viridis", alpha = 0.8) +
    scale_x_continuous(labels = scales::comma, breaks = scales::extended_breaks(n=7)) +
    labs(
      fill = ""
      , y = ""
      , x = "Fireshed Area (ha)"
    ) +
    theme_light() +
    theme(
      legend.position = "none" # c(0.9, 0.9)
      , axis.title = element_text(size=9)
      , axis.text.x = element_text(size=7)
    )

2.5.3 Map of Firesheds by Level of Constraint

# polish for mapping
map_fs_temp <- fireshed_proj_area %>% 
  dplyr::select(pa_id, geometry) %>% 
  dplyr::mutate(pa_id=as.character(pa_id)) %>% 
  dplyr::inner_join(
    constrained_by_scnro_ls_pa
    , by = dplyr::join_by("pa_id")
    , multiple = "all"
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ scales::percent(as.numeric(.x), accuracy = 0.1)
    )
    , dplyr::across(
      tidyselect::ends_with("_ha")
      , ~ scales::comma(as.numeric(.x), accuracy = 1)
    )
    , scenario_lab = scenario_lab %>% forcats::fct_rev()
  ) %>% 
  dplyr::mutate(
    Priority.Area.Name = area_name
    , Level.of.Constraint = cnstrnt_class
    , Pct.Treatable.Class = rmn_cnstrnt_class
    , Fireshed.ProjArea.ID = pa_id
    , Fireshed.ProjArea.Area.Hectares = feature_area_ha
    , Fireshed.WCS.Status = ifelse(is.na(fireshed_crisis_strategy),"Not High Risk",fireshed_crisis_strategy)
    , ForestShrub.Area.Hectares = covertype_area_ha
    , Pct.Rdctn.Protected = pct_rdctn1_protected
    , Pct.Rdctn.Slope = pct_rdctn2_slope
    , Pct.Rdctn.Roads = pct_rdctn3_roads
    , Pct.Rdctn.Riparian = pct_rdctn4_riparian
    , Pct.Rdctn.Administrative = pct_rdctn5_administrative
    , Treatable.Area.Hectares = rmn5_administrative_area_ha
    , Pct.Treatable.Area = pct_rmn5_administrative
  )
# filter for mapview
mapview_fs_temp <- map_fs_temp %>% 
  dplyr::filter(
    scenario_id == 2
  )
# ## filter for high priority only
# hi_map_fs_temp <- map_fs_temp %>% 
#   dplyr::filter(!is.na(fireshed_crisis_strategy))
# filter for mapview
mapview_fs_sc1_temp <- map_fs_temp %>% 
  dplyr::filter(
    scenario_id == 1
  )
mapview_fs_sc3_temp <- map_fs_temp %>% 
  dplyr::filter(
    scenario_id == 3
  )
mapview_fs_sc3_hi_temp <- map_fs_temp %>% 
  dplyr::filter(
    scenario_id == 3
    & fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
  )

# mapview
mapview::mapviewOptions(homebutton = FALSE, basemaps = c("OpenStreetMap","Esri.WorldImagery"))
mapview::mapview(wf_landscapes
        , color = "black"
        , lwd = 1
        , alpha.regions = 0
        , label = FALSE
        , legend = FALSE
        , popup = FALSE
        , layer.name = "WCS Landscapes"
        , hide = TRUE
  ) +
mapview::mapview(
  mapview_fs_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 2 Constraints"
  , hide = FALSE
  , popup = leafpop::popupTable(
      mapview_fs_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) + 
mapview::mapview(
  mapview_fs_sc1_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 1 Constraints"
  , hide = TRUE
  , popup = leafpop::popupTable(
      mapview_fs_sc1_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) + 
mapview::mapview(
  mapview_fs_sc3_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 3 Constraints"
  , hide = TRUE
  , popup = leafpop::popupTable(
      mapview_fs_sc3_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) +
mapview::mapview(
  mapview_fs_sc3_hi_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 3 Constraints (High Only)"
  , hide = TRUE
  , popup = leafpop::popupTable(
      mapview_fs_sc3_hi_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) 

# **Scenario 2 used for map above**

2.5.4 Map of Fireshed Project Areas by Scenario

plt_fn_temp <- function(scnum=1) {
  plt_basemap_simple +
    geom_sf(data = map_fs_temp %>% dplyr::filter(scenario_id==scnum)
            , aes(fill = cnstrnt_class), lwd = 0.05) +
    
  geom_sf(
    data = wf_landscapes
    # , mapping = aes(color = paste0(investment, " investment"))
    , fill=NA
    , color="black"
    , lwd=0.2
  )+
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = stringr::str_wrap(
          dplyr::case_when(
            stringr::str_count(name, "\\w+")<2
              ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
            , stringr::str_starts(name,"Sierra and Elko") ~ stringr::word(name,1,3)
            , TRUE ~ stringr::word(name,1,2)
          )
          , 7)
        # , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , color = "black"
    , stat = "sf_coordinates"
    , size = 2
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  # scale_color_manual(values = c("black", "gray33")) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  labs(
    fill = ""
    , tag = paste0(
      map_fs_temp %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::pull(scenario_lab) %>% unique()
      , " Map"
    )
  ) +
  theme(
    legend.position = c(0.8,0.9) # "top"
    , legend.text = element_text(size = 7)
    , plot.title = element_text(face = "bold", size = 9, margin = margin(0,0,30,0))
    , plot.tag = element_text(face = "bold", size = 10)
    , plot.tag.position = c(0.1,0.98) # "top"
    # , plot.tag.position = c(0.8,0.94) # "top"
  ) + 
  guides(color = "none")
}
# plot
# plt_fn_temp()
plt_fpas_sc_map <- 
  map_fs_temp %>% dplyr::pull(scenario_id) %>% unique() %>% 
    purrr::map(plt_fn_temp)
plt_fpas_sc_map

2.5.5 Map Fireshed Project Areas by Scenario and Landscape

How do small multiple maps of scenarios look?

the answer is: bad see this post about small multiple maps

area_nm_list <- sort(unique(constrained_by_scnro_ls$area_name))
facet_map_fn <- function(row_n=1) {
  plt <- ggplot() +
    geom_sf(
      data = map_fs_temp %>% 
        dplyr::filter(
          area_name == area_nm_list[row_n]
        )
      , aes(fill = cnstrnt_class), lwd = 0.05
    ) +
    geom_sf(
      data = wf_landscapes %>% 
        dplyr::filter(
          area_name == area_nm_list[row_n]
        )
      , color = "black"
      , fill=NA
      , lwd=0.2
    ) +
    facet_grid(
      cols = vars(scenario_lab)
      , rows = vars(area_name)
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
      , switch = "both"
      # , switch = "x"
    ) +
    scale_color_manual(values = c("black", "gray33")) +
    scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
    labs(
      fill = ""
      , tag = stringr::str_wrap(area_nm_list[row_n], 35)
        # scales::label_wrap(area_nm_list[row_n], 10)
    ) +
    theme_void() +
    # theme_light() +
    theme(
      legend.position = "none"
      # , legend.text = element_text(size = 7)
      , plot.subtitle = element_text(face = "bold", size = 9, hjust = 0.5, vjust = -5)
      # attempting to fix facet labels not aligning
      , strip.text.x = element_text(size = 8, vjust = 1)
      , strip.text.y = element_blank()
      , plot.tag.position = c(.5, 0.98)
      , plot.tag = element_text(size = 8, face = "bold")
      # , plot.background = element_rect(color = "gray88", fill=NA, size=0.5)
      , plot.margin = margin(0, 0, 0, 0, "cm")
      # facet labels not alinging in plot_grid
      # , strip.text.x = element_text(color = "black", face = "bold")
      # , strip.text.y.left = element_text(angle = 0)
      # , strip.placement = "outside"
      
    ) + 
    guides(color = "none")
  return(plt)
}
# print(facet_map_fn())
plt_list <- 1:length(area_nm_list) %>% 
  purrr::map(facet_map_fn)
# combine
scenario_maps <- cowplot::plot_grid(
  plotlist = plt_list #[1:4]
  , ncol = 3
  , align = "hv"
)
# scenario_maps
# plot_grid didn't work ;(
plt_list[16]
# plt_list

Make each plot individually by landscape

#########################################################################
#########################################################################
#########################################################################
##################hack to align plots for ggmap
ggmap_bbox_fn <- function(map, my_crs=3857) {
    if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
    # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
    # and set the names to what sf::st_bbox expects:
    map_bbox <- setNames(unlist(attr(map, "bb")), c("ymin", "xmin", "ymax", "xmax"))
    # Convert the bbox to an sf polygon, transform it to 3857, 
    # and convert back to a bbox (convoluted, but it works)
    bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), my_crs))
    # Overwrite the bbox of the ggmap object with the transformed coordinates 
    attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
    attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
    attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
    attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
    map
}
plt_crs <- 3857
#########################################################################
#########################################################################
#########################################################################
area_nm_list <- sort(unique(constrained_by_scnro_ls$area_name))
# buffer for smaller areas
if_smaller_ha_temp <- wf_landscapes %>% 
  dplyr::pull(hectares) %>% 
  quantile(0.45) %>% 
  purrr::pluck(1)
# create color palette
cols_rdylbu_lt <- c(
  "low constraint" = rev(RColorBrewer::brewer.pal(n=3,name="RdYlBu"))[1]
  , "med. constraint" = rev(RColorBrewer::brewer.pal(n=3,name="RdYlBu"))[2]
  , "high constraint" = rev(RColorBrewer::brewer.pal(n=3,name="RdYlBu"))[3]
)
# function to plot basemap
landscape_basemap_fn <- function(row_n=1) {
  # should zoom in?
  zoom_level <- dplyr::case_when(
    area_nm_list[row_n] %in% c(
      "WA: Colville Northeast Washington Vision"
      , "CA: Southern California Fireshed Risk Reduction Strategy"
      ) ~ 7
    , area_nm_list[row_n] %in% c(
      "UT: Wasatch"
      , "CA: Trinity Forest Health and Fire-Resilient Rural Communities"
      , "ID: Nez Perce-Clearwater-Lower Salmon"
      , "NV: Sierra and Elko Fronts"
      , "OR: Mount Hood Forest Health and Fire-Resilient Communities"
    ) ~ 8
    , area_nm_list[row_n] %in% c(
      "CA: Plumas Community Protection"
      , "CA: North Yuba"
      , "MT: Kootenai Complex"
      , "UT: Pine Valley"
    ) ~ 9
    , 
      (wf_landscapes %>% 
        dplyr::filter(
            area_name == area_nm_list[row_n]
        ) %>% 
        dplyr::pull(hectares) %>% 
        purrr::pluck(1)
        <= if_smaller_ha_temp
      ) ~ 10
    , TRUE ~ 8
  )
  # should buffer extend?
  buffer_box <- dplyr::case_when(
    area_nm_list[row_n] %in% c("WA: Colville Northeast Washington Vision") ~ 50000
    , area_nm_list[row_n] %in% c(
      "CA: Southern California Fireshed Risk Reduction Strategy"
      , "ID: Nez Perce-Clearwater-Lower Salmon"
      , "MT: Kootenai Complex"
    ) ~ 40000
    , area_nm_list[row_n] %in% c(
      "AZ: San Carlos Apache Tribal Forest Protection"
      , "NV: Sierra and Elko Fronts"
      , "UT: Pine Valley"
    ) ~ 21000
    , area_nm_list[row_n] %in% c(
      "CA: Plumas Community Protection"
      , "CA: North Yuba"
      , "CA: Trinity Forest Health and Fire-Resilient Rural Communities"
      , "OR: Mount Hood Forest Health and Fire-Resilient Communities"
    ) ~ 35000
    , area_nm_list[row_n] %in% c("UT: Wasatch") ~ 5000
    , 
      (wf_landscapes %>% 
        dplyr::filter(
            area_name == area_nm_list[row_n]
        ) %>% 
        dplyr::pull(hectares) %>% 
        purrr::pluck(1)
        <= if_smaller_ha_temp
      ) ~ 15000
    , TRUE ~ 5000
  )
  
  # bounding box
  bb_temp <- 
    map_fs_temp %>% 
    dplyr::filter(
      area_name == area_nm_list[row_n]
    ) %>% 
    sf::st_union() %>% 
    sf::st_transform(crs=5070) %>% 
    sf::st_buffer(as.numeric(buffer_box)) %>% 
    sf::st_transform(crs=4326) %>% # same as get_map return
    sf::st_bbox()
  # set bbox for get call
  bbox_temp <- c(
    bottom = bb_temp[[2]]
    , top = bb_temp[[4]]
    , right = bb_temp[[3]]
    , left = bb_temp[[1]]
  )
  # get map
  hey_ggmap <- ggmap::get_stamenmap(
    bbox = bbox_temp
    , zoom = zoom_level
    , maptype = "terrain" #"toner-hybrid" # "terrain"
    , crop = T
  )
  # ggmap(hey_ggmap)
  # apply align function
    hey_ggmap_aligned <- ggmap_bbox_fn(hey_ggmap, plt_crs) # Use the function
  # plot
  plt <- ggmap(hey_ggmap_aligned) + 
    coord_sf(
      expand = FALSE
    ) +
    labs(
      title = area_nm_list[row_n]
    ) +
    theme_light() +
    theme(
      legend.position = "bottom"
      , legend.direction  = "horizontal"
      , legend.title = element_text(size=7)
      , legend.margin = margin(c(-10,0,0,0))
      , plot.title = element_text(size = 11)
      , strip.text = element_text(color = "black", face = "bold")
      , axis.title = element_blank()
      , axis.text = element_blank()
      , axis.ticks = element_blank()
      , panel.grid = element_blank()
      , plot.margin = margin(0, 0, 0, 0, "cm")
    )
  # return
  return(plt)
}
# landscape_basemap_fn(row_n = 9)
# function to plot scenario constraints + basemap
scenario_cnstrnt_map_fn <- function(row_n = 1) {
  plt <- landscape_basemap_fn(row_n) + 
      geom_sf(
        data = map_fs_temp %>%
          dplyr::filter(
            area_name == area_nm_list[row_n]
          ) %>%
          sf::st_transform(crs=plt_crs)
        , aes(fill = cnstrnt_class), lwd = 0.05
        , inherit.aes = F
        , alpha = 0.6
      ) +
      geom_sf(
        data = wf_landscapes %>% 
          dplyr::filter(
              area_name == area_nm_list[row_n]
            ) %>% 
          sf::st_transform(crs=plt_crs)
        , fill = NA, color = "black", lwd = 0.3
        , inherit.aes = F
      ) +
      facet_grid(
        cols = vars(scenario_lab)
        , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
        # , switch = "both"
      ) +
      scale_fill_manual(values = cols_rdylbu_lt) +
      labs(
        fill = "Level of\nConstraint"
      ) +
        guides(
        fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
      )
  # return
  return(plt)
}
# scenario_cnstrnt_map_fn(row_n = 9)

# go through list of areas to make plot
plt_list <- 1:length(area_nm_list) %>%
  purrr::map(scenario_cnstrnt_map_fn)
# scenario_maps
plt_list[]

####
# area_nm_list
# scenario_cnstrnt_map_fn(18)

2.5.6 Patch Analysis

Notes from Collins et al. 2010:

The larger the individual fuel treatment the greater the potential for tree survival, because of reduced edge effect mortality (Ritchie et al. 2007). Larger individual treatments have a greater potential to reduce fire behavior and slow fire spread, which ultimately impacts adjacent untreated stands and should enhance suppression opportunities and increase firefighter safety (Collins et al. 2010). Planning fewer, larger individual treatments across the landscape appears to be a better strategy when human community protection is a primary concern (Schmidt et al. 2008). Finney et al. (2007) note that the relative improvement of optimized treatment placement breaks down when larger proportions of the landscape (40–50%) are excluded from treatment because of land-management constraints.

# create spatial data of interconnected patches by constraint class
# this is to represent the spatial aggregation of fireshed project areas of similar levels of  mechanical constraint
  # aggregation = Tendency of patch or land-cover types to be spatially adjacent or in close proximity
cnstrnt_class_patches <-
  fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::select(pa_id, area_name, scenario_id, scenario_lab, cnstrnt_class)
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      ) %>% 
  # union vectors/polygons by constraint class into one multipolygon
    # this dissolves inner boundaries by constraint class
  dplyr::group_by(area_name, scenario_id, scenario_lab, cnstrnt_class) %>% 
  dplyr::summarize(
    geometry = sf::st_union(geometry, by_feature = F)
  ) %>% 
  # separate a multipolygon geometry into several polygons objects after performing a st_union()
  sf::st_cast(to = "MULTIPOLYGON") %>%
  st_cast(to = "POLYGON") %>% 
  # calculate patch area 
  dplyr::ungroup() %>% 
  dplyr::group_by(area_name, scenario_id, scenario_lab, cnstrnt_class) %>% 
  dplyr::mutate(
    patch_area_ha = as.numeric(sf::st_area(geometry))/10000
    , patch_area_rank = dplyr::dense_rank(dplyr::desc(patch_area_ha))
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::arrange(scenario_id, area_name, cnstrnt_class, dplyr::desc(patch_area_ha)) %>% 
  # calculate proportion of total area in each patch
  dplyr::inner_join(
    # total area of fireshed project areas included for analysis
      # after union to remove internal boundaries
    fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(scenario_id==1) %>% 
          dplyr::select(pa_id, area_name)
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      ) %>% 
      dplyr::group_by(area_name) %>% 
      dplyr::summarize(
        geometry = sf::st_union(geometry, by_feature = F)
      ) %>%
      dplyr::mutate(
        total_area_ha = as.numeric(sf::st_area(geometry))/10000
      ) %>% 
      sf::st_drop_geometry()
    , by = dplyr::join_by(area_name)
  ) %>% 
  dplyr::mutate(
    pct_patch_landscape = patch_area_ha / total_area_ha
  )

# function to highlight largest patch with fireshed project areas under
cols_rdylbu_dk <- c(
  "low constraint" = rev(RColorBrewer::brewer.pal(n=11,name="RdYlBu"))[2]
  , "med. constraint" = rev(RColorBrewer::brewer.pal(n=11,name="RdYlBu"))[7]
  , "high constraint" = rev(RColorBrewer::brewer.pal(n=11,name="RdYlBu"))[10]
)
largest_patch_map_fn <- function(nrow = 1){
    plt <- scenario_cnstrnt_map_fn(nrow) +
      geom_sf(
          data = cnstrnt_class_patches %>%
            dplyr::filter(
              area_name == area_nm_list[nrow]
              & patch_area_rank == 1
              & cnstrnt_class != "med. constraint"
            ) %>%
            sf::st_transform(crs=plt_crs)
          , aes(color = cnstrnt_class, fill = cnstrnt_class), lwd = 1.4
          , inherit.aes = F
          , alpha = 0
          # , show.legend = F
        ) +
      scale_color_manual(values = cols_rdylbu_dk, drop = F) +
      labs(
        color = "Level of\nConstraint"
        , subtitle = "largest patch of contiguous fireshed project areas classified by level of mechanical constraint" # largest patch of spatially adjacent...
      ) +
      theme(
        plot.subtitle = element_text(size = 8)
      ) +
      guides(
        color = "none"
        , fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9, linewidth = 0))
      )
  # return
  return(plt)
}
# largest_patch_map_fn(13)

# plot each landscape
plt_ptch_list <- 1:length(area_nm_list) %>%
  purrr::map(largest_patch_map_fn)
# scenario_maps
plt_ptch_list[]

Comparison of largest patch of lightly constrained fireshed project areas versus largest patch of highly constrained fireshed project areas.

Planning fewer, larger individual treatments across the landscape appears to be a better strategy when human community protection is a primary concern (Schmidt et al. 2008). Finney et al. (2007) note that the relative improvement of optimized treatment placement breaks down when larger proportions of the landscape (40–50%) are excluded from treatment because of land-management constraints.

# aggregate to constraint class level
cnstrnt_class_patches %>% 
  sf::st_drop_geometry() %>% 
  dplyr::group_by(area_name, scenario_lab, cnstrnt_class) %>% 
  dplyr::summarize(
    max_patch_area_ha = max(patch_area_ha, na.rm = T)
    , mean_patch_area_ha = mean(patch_area_ha, na.rm = T)
    , max_pct_patch_landscape = max(pct_patch_landscape, na.rm = T)
    , sum_pct_patch_landscape = sum(pct_patch_landscape)
  ) %>% 
  dplyr::ungroup() %>%
  dplyr::mutate(
    max_pct_patch_landscape_hack = dplyr::case_when(
      cnstrnt_class=="high constraint" ~ -max_pct_patch_landscape
      , TRUE ~ max_pct_patch_landscape
    )
  ) %>% 
  dplyr::filter(
    cnstrnt_class!="med. constraint"
  ) %>% 
# plot
ggplot(
  mapping = aes(
    y = scenario_lab
    , x = max_pct_patch_landscape_hack
    , fill = cnstrnt_class
    , group = cnstrnt_class
  )
) +
  # geom_line(
  #   mapping = aes(x = forcats::fct_rev(scenario_lab), y = max_pct_patch_landscape, color = forcats::fct_rev(cnstrnt_class), group = forcats::fct_rev(cnstrnt_class))
  #   , linewidth = 1.2
  # ) +
  geom_col(
    width = 0.7, alpha = 0.9
  ) +
  geom_vline(xintercept = 0, color = "gray55", linetype = "solid", linewidth = 1.1) +
  geom_text(
    mapping = aes(
      label = ifelse(
        max_pct_patch_landscape_hack<0
        , scales::percent(max_pct_patch_landscape, accuracy = 1)
        , ""
      )
      , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
      , fontface = "bold"
    )
    , color = "black", size = 2.7
    , hjust = +1
  ) +
  geom_text(
    mapping = aes(
      label = ifelse(
        max_pct_patch_landscape_hack>=0
        , scales::percent(max_pct_patch_landscape, accuracy = 1)
        , ""
      )
      , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
      , fontface = "bold"
    )
    , color = "black", size = 2.7
    , hjust = 0
  ) +
  annotate(
    geom = "text"
    , x = -0.77
    , y = 0
    , label = "high constraint"
    , color = cols_rdylbu_lt[3]
    , fontface = "bold"
    , size = 2.5
  ) +
  annotate(
    geom = "text"
    , x = 0.86
    , y = 0
    , label = "low constraint"
    , color = cols_rdylbu_lt[1]
    , fontface = "bold"
    , size = 2.5
  ) +
  scale_fill_manual(values = cols_rdylbu_lt) +
  scale_x_continuous(
    limits = c(-1.05,1.13)
    , labels = scales::percent_format()
    # , expand = expansion(mult = c(0, 0.1))
    # , breaks = scales::extended_breaks(5)
  ) +
  scale_y_discrete(expand = expansion(mult = c(0.8,0.3))) +
  facet_wrap(
    facets = vars(area_name)
    , ncol = 3
    , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    , scales = "free_y"
  ) +
  labs(
    fill = "Level of\nConstraint"
    , subtitle = "percentage of total landscape area comprised by the largest patch" # "Largest Patch % of Landscape Area"
    , x = "" # "Largest Patch % of Landscape Area"
    , y = ""
  ) +
  theme_light() +
  theme(
    legend.position = "none" # "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.ticks.x = element_blank()
    , panel.grid.major = element_blank()
    , panel.grid.minor = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 9)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

2.6 High-risk fireshed project areas

Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the plan (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).

2.6.1 Project Areas by WCS Status

# aggregate data by area, status
wcs_status_temp <- constrained_by_scnro_ls_pa %>% 
  dplyr::filter(scenario_id==1) %>%
  dplyr::group_by(area_name,fireshed_crisis_strategy) %>% 
  dplyr::summarise(
    n = n()
    , feature_area_ha = sum(feature_area_ha)
  ) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(
    tot_n = sum(n)
    , tot_area = sum(feature_area_ha)
    , pct_n = n/tot_n
    , pct_area = feature_area_ha/tot_area
  ) %>% 
  dplyr::ungroup()
# plot
ggplot(data = wcs_status_temp) +
  geom_col(
    mapping = aes(y = reorder(area_name,tot_n), x = n, fill = fireshed_crisis_strategy)
    , color = NA, width = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = area_name, x = n, group = fireshed_crisis_strategy
      , label = scales::comma(ifelse(n>20,n,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_pa %>% 
        dplyr::filter(scenario_id==1) %>%
        dplyr::group_by(area_name) %>% 
        dplyr::summarise(
          n = n()
          , feature_area_ha = sum(feature_area_ha)
        )
    , mapping = aes(
      y = reorder(area_name,n), x = n
      , label = scales::comma(n, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray")) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.1)),labels = scales::comma_format()) +
  labs(
    fill = "High-Risk\nStatus"
    , y = ""
    , x = "# of Fireshed\nProject Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
  )

2.6.2 Area of high-risk project areas

Only including fireshed project areas that have at least 25% of their area within the bounds of the priority landscapes and are in a high-risk fireshed.

See this section for discussion of high-risk firesheds and criteria for including fireshed project areas.

# data filter
wcs_status_temp %>% 
  dplyr::filter(fireshed_crisis_strategy!="Not High Risk") %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(tot_area = sum(feature_area_ha)) %>% 
# plot
ggplot() +
  geom_col(
    mapping = aes(y = reorder(area_name,tot_area), x = feature_area_ha, fill = fireshed_crisis_strategy)
    , color = NA, width = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = area_name, x = feature_area_ha, group = fireshed_crisis_strategy
      , label = ifelse(
          feature_area_ha<150000
          ,NA
          ,scales::comma(feature_area_ha,suffix = "k", scale = 1e-3, accuracy = 1)
        )
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_pa %>% 
        dplyr::filter(
          scenario_id==1
          & fireshed_crisis_strategy!="Not High Risk"
        ) %>%
        dplyr::group_by(area_name) %>% 
        dplyr::summarise(
          n = n()
          , feature_area_ha = sum(feature_area_ha)
        )
    , mapping = aes(
      y = reorder(area_name,feature_area_ha), x = feature_area_ha
      , label = scales::comma(feature_area_ha,suffix = "k ha", scale = 1e-3, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray")) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.15)),labels = scales::comma_format(suffix = "k ha", scale = 1e-3, accuracy = 1)) +
  labs(
    fill = "High-Risk\nStatus"
    , y = ""
    , x = "Area (ha) of Fireshed\nProject Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
  )

2.6.3 Reduction by Constraint

See this section for reduction by constraint at the landscape-level.

# aggregate by area_name, scenario for hi-risk pa's
constrained_by_scnro_ls_hionly_long <-
  constrained_by_scnro_ls_pa %>%
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::select(scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
                , (tidyselect::starts_with("rmn") & tidyselect::ends_with("area_ha"))
                , covertype_area_ha
  ) %>%
  dplyr::group_by(
    scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
  ) %>% 
  dplyr::summarise(
    dplyr::across(
      tidyselect::ends_with("area_ha")
      , sum
    )
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("rmn")
    , names_to = "constraint"
    , values_to = "rmn_area_ha"
    , names_prefix = "rmn"
    , values_drop_na = F
  ) %>%
  dplyr::mutate(constraint = stringr::str_remove(constraint,"_area_ha")) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , constraint = factor(
        constraint
        , ordered = TRUE
        , levels = c(
          "protected"
          , "slope"
          , "roads"
          , "riparian"
          , "administrative"
        )
        , labels = c(
          "Protected"
          , "Slope\nSteepness"
          , "Road\nDistance"
          , "Riparian\nBuffer"
          , "Administrative"
        )
      ) %>% forcats::fct_rev()
    , pct_rmn = rmn_area_ha/covertype_area_ha
    , pct_rdctn = dplyr::case_when(
      constraint_lvl == 1 ~ -1*(1-pct_rmn)
      , TRUE ~ -1*(dplyr::lag(pct_rmn)-pct_rmn)
    )
    , pct_rdctn_total = max(-1*(1 - ifelse(constraint_lvl == 5, pct_rmn,0)))
  )
  
# plot
ggplot(data = constrained_by_scnro_ls_hionly_long) +
  geom_col(
    mapping = aes(y = scenario_lab, x = pct_rdctn, fill = constraint)
    , color = NA, width = 0.7
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = scenario_lab, x = pct_rdctn
      , label = scales::percent(ifelse(pct_rdctn<.06*-1,pct_rdctn,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_hionly_long %>% dplyr::filter(constraint_lvl==5)
    , mapping = aes(
      y = scenario_lab, x = -1*(1 - pct_rmn)
      , label = scales::percent(-1*(1 - pct_rmn), accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_wrap(facets = vars(area_name)
      , ncol = 3
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    ) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0.02, .21)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
    , subtitle = "High-Risk Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 10) 
  )

The figure above shows the reduction by constraint for only high-risk fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape.

2.6.4 Reduction by Constraint high-risk vs. total

# combine data and plot
pct_rdctn_temp <- 
constrained_by_scnro_ls_hionly_long %>% 
  dplyr::ungroup() %>% 
  dplyr::select(scenario_lab, area_name, pct_rdctn, constraint_lvl) %>% 
  dplyr::mutate(dta_src = "High-Risk Project Areas") %>% 
  dplyr::bind_rows(
    constrained_by_scnro_ls_long %>% 
      dplyr::ungroup() %>% 
      dplyr::select(scenario_lab, area_name, pct_rdctn, constraint_lvl) %>% 
      dplyr::filter(!is.na(constraint_lvl)) %>% 
      dplyr::mutate(dta_src = "Overall Landscape Area")
  ) %>% 
  dplyr::group_by(scenario_lab, area_name, dta_src) %>% 
  dplyr::mutate(
    pct_rdctn_total = sum(pct_rdctn)
    , constraint = factor(
        constraint_lvl
        , ordered = TRUE
        , levels = 1:5
        , labels = c(
          "Protected"
          , "Slope\nSteepness"
          , "Road\nDistance"
          , "Riparian\nBuffer"
          , "Administrative"
        )
      ) %>% forcats::fct_rev()
    , scenario_lab = scenario_lab %>% forcats::fct_rev()
  )
# plot
ggplot(data = pct_rdctn_temp) +
  geom_col(
    mapping = aes(y = dta_src, x = pct_rdctn, fill = constraint)
    , color = NA, width = 0.7
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = dta_src, x = pct_rdctn
      , label = scales::percent(ifelse(pct_rdctn<.075*-1,pct_rdctn,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  geom_text(
    data = pct_rdctn_temp %>% dplyr::filter(constraint_lvl==5)
    , mapping = aes(
      y = dta_src, x = pct_rdctn_total
      , label = scales::percent(pct_rdctn_total, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_grid(
      cols = vars(scenario_lab)
      , rows = vars(area_name)
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
      , switch = "y"
    ) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0.02, .21)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_blank()
    , strip.text.x = element_text(color = "black", face = "bold")
    , strip.text.y.left = element_text(angle = 0)
    , strip.placement = "outside"
    
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

2.7 Objective vs. Treatable

Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the Strategy (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).

Using the area of fireshed project areas that have at least 25% of their area within the bounds of the priority landscapes and are in a high-risk fireshed as the treatment area objective (20-40%).

# Area of pa's by level of constraint
constrained_by_scnro_ls_cnstrnt_temp <- constrained_by_scnro_ls_pa %>% 
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::group_by(scenario_id,scenario_lab,state,name,area_name,cnstrnt_class) %>% 
  dplyr::summarise(
    n = n()
    , dplyr::across(
      c(feature_area_ha,covertype_area_ha,rmn5_administrative_area_ha)
      , ~ sum(.x,na.rm = T)
    )
  ) %>% 
  dplyr::group_by(scenario_id,scenario_lab,state,name,area_name) %>% 
  dplyr::mutate(
    # cnstrnt_class = stringr::word(cnstrnt_class) %>% stringr::str_remove_all("\\.")
    total_area_ha = sum(feature_area_ha)
    , total_covertype_area_ha = sum(covertype_area_ha)
    # "the need to treat 20 to 40 percent of the overall fireshed"
    , pct_treatable =  rmn5_administrative_area_ha/total_area_ha
    , objective_lo_treat_area_ha = total_area_ha*.2
    , objective_hi_treat_area_ha = total_area_ha*.4
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::arrange(state,name,scenario_id,cnstrnt_class)

2.7.1 Objective vs. Treatable Table

table_temp <-
  constrained_by_scnro_ls_cnstrnt_temp %>% 
  dplyr::mutate(
    cnstrnt_class = stringr::word(cnstrnt_class) %>% stringr::str_remove_all("\\.")
  ) %>% 
  dplyr::rename(
    remaining_pct=pct_treatable
    , remaining_ha=rmn5_administrative_area_ha
  ) %>%
  dplyr::select(-c(feature_area_ha, covertype_area_ha)) %>% 
  tidyr::pivot_wider(
    names_from = cnstrnt_class
    , values_from = c(n, remaining_pct, remaining_ha)
    , names_glue = "{cnstrnt_class}_{.value}"
    , values_fill = 0
  ) %>% 
  dplyr::mutate(
    tot_remaining_ha = low_remaining_ha+med_remaining_ha
    ,tot_remaining_pct = low_remaining_pct+med_remaining_pct
  ) %>% 
  dplyr::relocate(tot_remaining_ha,.after = med_remaining_ha) %>% 
  dplyr::inner_join(
    wf_landscapes %>% 
      sf::st_drop_geometry() %>% 
      dplyr::select(area_name,hectares) %>% 
      dplyr::rename(landscape_area_ha=hectares)
    , by = dplyr::join_by(area_name)
  )
# make table 
kableExtra::kable(
    table_temp %>% 
      dplyr::select(
        c(
          scenario_lab
          # , landscape_area_ha
          , total_area_ha
          , total_covertype_area_ha
          # light and mod area
          , low_remaining_ha, med_remaining_ha, tot_remaining_ha
          , low_remaining_pct, med_remaining_pct, tot_remaining_pct
        )
      ) %>% 
      dplyr::mutate(
        dplyr::across(
          tidyselect::ends_with("_ha")
          , ~ scales::comma(.x, suffix = "k", scale = 1e-3, accuracy = 1)
        )
        , dplyr::across(
          tidyselect::ends_with("_pct")
          , ~ scales::percent(.x, accuracy = .1)
        )
      )
  , caption = "Wildfire Crisis Strategy High-Risk Fireshed Project Areas<br>mechanically available for treatment in lightly and moderately constrained project areas (≥25% in landscape)"
  , col.names = c(
    ""
    # , "Total\nLandscape"
    , "Total Area (ha)"
    , "Forest+Shrub Area (ha)"
    , "Lightly Constrained", "Moderately Constrained"
    , "Total"
    , "Lightly Constrained", "Moderately Constrained"
    , "Total"
  
  )
  , escape = F
) %>%  
    add_header_above(c(" " = 1, "High-Risk Fireshed\nProject Areas" = 2, "Treatable High-Risk\nArea (ha)"=3, "Treatable High-Risk\nPercent (%)" = 3)) %>%
    kable_classic(full_width=T) %>% 
    pack_rows(index = table(forcats::fct_inorder(table_temp$area_name))) %>% 
    kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%  
    kableExtra::column_spec(
      9 # 10
      , color = "white"
      , background = 
          table_temp %>% 
            dplyr::mutate(
              is_gt20 = dplyr::case_when(
                low_remaining_pct >= .2 ~ "#008b8b"
                , tot_remaining_pct >= .2 ~ "#545454"
                , TRUE ~ "#cd3700"
              )
          ) %>% 
          dplyr::pull(is_gt20)
          
    ) %>% 
    kableExtra::scroll_box(width = "740px")

Table 2.3: Wildfire Crisis Strategy High-Risk Fireshed Project Areas
mechanically available for treatment in lightly and moderately constrained project areas (≥25% in landscape)

	High-Risk Fireshed Project Areas		Treatable High-Risk Area (ha)			Treatable High-Risk Percent (%)
	Total Area (ha)	Forest+Shrub Area (ha)	Lightly Constrained	Moderately Constrained	Total	Lightly Constrained	Moderately Constrained	Total
AZ: 4FRI
Scenario 1	1,579k	1,436k	360k	85k	445k	22.8%	5.4%	28.2%
Scenario 2	1,579k	1,436k	554k	58k	611k	35.0%	3.7%	38.7%
Scenario 3	1,579k	1,436k	936k	30k	966k	59.3%	1.9%	61.2%
AZ: Prescott
Scenario 1	221k	191k	34k	19k	53k	15.3%	8.6%	23.9%
Scenario 2	221k	191k	69k	14k	83k	31.3%	6.2%	37.6%
Scenario 3	221k	191k	112k	0k	112k	50.9%	0.0%	50.9%
AZ: San Carlos Apache Tribal Forest Protection
Scenario 1	250k	230k	4k	29k	33k	1.6%	11.6%	13.2%
Scenario 2	250k	230k	42k	36k	78k	16.8%	14.4%	31.2%
Scenario 3	250k	230k	80k	20k	100k	31.9%	7.9%	39.8%
CA: North Yuba
Scenario 1	206k	197k	32k	26k	59k	15.8%	12.8%	28.6%
Scenario 2	206k	197k	85k	6k	91k	41.4%	2.8%	44.2%
Scenario 3	206k	197k	104k	7k	111k	50.5%	3.6%	54.1%
CA: Plumas Community Protection
Scenario 1	157k	148k	24k	25k	49k	15.3%	15.7%	31.0%
Scenario 2	157k	148k	68k	5k	72k	43.1%	3.0%	46.1%
Scenario 3	157k	148k	91k	0k	91k	57.9%	0.0%	57.9%
CA: Southern California Fireshed Risk Reduction Strategy
Scenario 1	1,053k	854k	44k	80k	124k	4.2%	7.6%	11.8%
Scenario 2	1,053k	854k	186k	68k	255k	17.7%	6.5%	24.2%
Scenario 3	1,053k	854k	253k	65k	318k	24.0%	6.1%	30.2%
CA: Stanislaus
Scenario 1	153k	143k	31k	21k	52k	20.4%	13.7%	34.1%
Scenario 2	153k	143k	63k	6k	70k	41.4%	4.1%	45.5%
Scenario 3	153k	143k	84k	0k	84k	54.9%	0.0%	54.9%
CA: Trinity Forest Health and Fire-Resilient Rural Communities
Scenario 1	394k	336k	0k	30k	30k	0.0%	7.5%	7.5%
Scenario 2	394k	336k	32k	37k	68k	8.1%	9.3%	17.4%
Scenario 3	394k	336k	140k	20k	160k	35.6%	5.0%	40.6%
CO: Colorado Front Range
Scenario 1	1,138k	823k	120k	100k	220k	10.5%	8.8%	19.3%
Scenario 2	1,138k	823k	301k	57k	358k	26.5%	5.0%	31.4%
Scenario 3	1,138k	823k	398k	43k	441k	34.9%	3.8%	38.7%
ID: Nez Perce-Clearwater-Lower Salmon
Scenario 1	520k	444k	91k	39k	130k	17.4%	7.6%	25.0%
Scenario 2	520k	444k	203k	12k	215k	39.1%	2.2%	41.3%
Scenario 3	520k	444k	213k	8k	221k	40.9%	1.6%	42.5%
ID: Southwest Idaho
Scenario 1	734k	632k	243k	32k	275k	33.1%	4.3%	37.4%
Scenario 2	734k	632k	408k	20k	428k	55.6%	2.7%	58.3%
Scenario 3	734k	632k	429k	11k	441k	58.4%	1.5%	60.0%
MT: Kootenai Complex
Scenario 1	322k	294k	101k	13k	114k	31.5%	4.0%	35.5%
Scenario 2	322k	294k	143k	6k	149k	44.5%	1.9%	46.4%
Scenario 3	322k	294k	194k	7k	200k	60.2%	2.0%	62.3%
NV: Sierra and Elko Fronts
Scenario 1	678k	453k	119k	39k	157k	17.5%	5.7%	23.2%
Scenario 2	678k	453k	231k	17k	248k	34.1%	2.5%	36.7%
Scenario 3	678k	453k	257k	13k	270k	37.9%	2.0%	39.9%
NM: Enchanted Circle
Scenario 1	318k	285k	51k	29k	79k	16.0%	9.0%	25.0%
Scenario 2	318k	285k	122k	13k	135k	38.3%	4.1%	42.5%
Scenario 3	318k	285k	142k	6k	149k	44.7%	2.0%	46.7%
OR: Central Oregon
Scenario 1	278k	232k	145k	5k	150k	52.2%	1.8%	54.0%
Scenario 2	278k	232k	169k	2k	170k	60.5%	0.6%	61.2%
Scenario 3	278k	232k	203k	0k	203k	73.1%	0.0%	73.1%
OR: Klamath River Basin
Scenario 1	971k	808k	33k	63k	96k	3.4%	6.5%	9.8%
Scenario 2	971k	808k	190k	47k	237k	19.5%	4.9%	24.4%
Scenario 3	971k	808k	316k	48k	364k	32.5%	4.9%	37.4%
OR: Mount Hood Forest Health and Fire-Resilient Communities
Scenario 1	339k	244k	22k	22k	44k	6.4%	6.6%	13.0%
Scenario 2	339k	244k	49k	21k	71k	14.6%	6.3%	20.9%
Scenario 3	339k	244k	135k	8k	143k	39.8%	2.3%	42.1%
UT: Pine Valley
Scenario 1	187k	173k	18k	20k	39k	9.9%	10.9%	20.8%
Scenario 2	187k	173k	53k	12k	65k	28.5%	6.3%	34.9%
Scenario 3	187k	173k	57k	13k	69k	30.3%	6.8%	37.2%
UT: Wasatch
Scenario 1	156k	126k	23k	10k	33k	14.8%	6.6%	21.4%
Scenario 2	156k	126k	45k	8k	52k	28.6%	5.0%	33.7%
Scenario 3	156k	126k	50k	6k	56k	31.9%	4.0%	36.0%
WA: Central Washington Initiative
Scenario 1	818k	588k	21k	32k	52k	2.5%	3.9%	6.4%
Scenario 2	818k	588k	76k	44k	120k	9.3%	5.3%	14.6%
Scenario 3	818k	588k	258k	29k	287k	31.6%	3.6%	35.1%
WA: Colville Northeast Washington Vision
Scenario 1	272k	244k	70k	13k	83k	25.9%	4.7%	30.6%
Scenario 2	272k	244k	122k	8k	129k	44.7%	2.9%	47.6%
Scenario 3	272k	244k	149k	6k	155k	54.7%	2.2%	56.9%

2.7.2 Reduction in treatable area flow

table_temp %>% 
  dplyr::filter(scenario_id==1) %>% 
  dplyr::select(
    c(
      area_name
      , landscape_area_ha
      , total_area_ha
      , total_covertype_area_ha
      # , tot_remaining_ha
      # , low_remaining_ha
    )
  ) %>%
  tidyr::pivot_longer(
    !area_name
    , names_to = "name"
    , values_to = "hectares"
  ) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(order = dplyr::row_number(),y_val=2) %>% 
  dplyr::ungroup() %>% 
  dplyr::cross_join(data.frame(scenario_id = 1:3)) %>% 
  ####
  dplyr::bind_rows(
    table_temp %>% 
      dplyr::select(
        c(
          area_name
          , scenario_id
          , tot_remaining_ha
          , low_remaining_ha
        )
      ) %>%
      tidyr::pivot_longer(
        !c(area_name,scenario_id)
        , names_to = "name"
        , values_to = "hectares"
      ) %>% 
      dplyr::group_by(area_name,scenario_id) %>% 
      dplyr::mutate(
        order = dplyr::row_number()+3
        , y_val = scenario_id
      ) %>% 
      dplyr::ungroup() %>% 
      dplyr::relocate(scenario_id,.after = dplyr::last_col())
  ) %>% 
  dplyr::arrange(area_name,order,y_val ) %>% 
  dplyr::mutate(
    order_lab = factor(
      order
      , levels = 1:5
      , labels = c(
        "Total\nLandscape"
        , "High-risk\nFireshed PAs"
        , "High-risk\nForest+Shrub"
        , "Lt.+Md. Constrained"
        , "Lightly Constrained"
      )
      , ordered = T
    )
  ) %>% 
  dplyr::inner_join(
    wf_landscapes %>% 
      sf::st_drop_geometry() %>% 
      dplyr::select(area_name, hectares) %>% 
      dplyr::rename(tot_hectares=hectares)
    , by = dplyr::join_by(area_name)
  ) %>% 
  dplyr::inner_join(
    table_temp %>% 
      dplyr::filter(scenario_id==1) %>% 
      dplyr::select(area_name, objective_lo_treat_area_ha)
    , by = dplyr::join_by(area_name)
  ) %>% 
  dplyr::mutate(
    y_val = -y_val
    , sizing = ifelse(hectares>tot_hectares,1,hectares/tot_hectares)
    , coloring = dplyr::case_when(
        order %in% 4:5
          & hectares >= objective_lo_treat_area_ha
          ~ 1
        , order %in% 4:5
          & hectares < objective_lo_treat_area_ha
          ~ 2
        , T ~ 3
    )
  ) %>% 
  # dplyr::filter(stringr::str_starts(area_name,"CA")) %>% 
  # View()
ggplot() +
  geom_line(
    mapping = aes(x = order_lab, y = y_val, group = scenario_id)
    , color = "gray"
  ) +
  geom_point(
    mapping = aes(x = order_lab, y = y_val, size = sizing*100, color = factor(coloring))
    , shape = 15
  ) +
  geom_text(
    mapping = aes(
      x = order_lab, y = y_val, group = scenario_id
      , label = scales::comma(hectares,suffix = "k", scale = 1e-3, accuracy = 1)
    )
    , size = 2.5
    , color = "black"
    , vjust = 2.5
    , angle = 90
  ) +
  geom_text(
    mapping = aes(
      x = order_lab, y = y_val, group = scenario_id
      , label = ifelse(
        order == 4, paste0("Scenario", y_val*-1), NA
      )
      , fontface = "italic"
    )
    , size = 2.5
    , color = "black"
    , hjust = 0.1
    , vjust = -1.4
  ) +
  facet_wrap(facets = vars(area_name)
    , ncol = 3
    , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    , scales = "free_x"
  ) +
  scale_color_manual(values = c("#008b8b", "#cd3700", "#bebebe")) +
  scale_y_continuous(expand = expansion(c(0.3,0.3))) +
  scale_x_discrete(labels = scales::label_wrap(10)) +
  labs(
    x = ""
    , y = ""
    , subtitle = "<span><span style='color:#008b8b;'><b><i>≥20% treatable</i></b></span> | <span style='color:#cd3700;'><b><i><20% treatable</i></b></span></span>"
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.y = element_blank()
    , axis.ticks.y = element_blank()
    , axis.text = element_text(size = 7)
    , panel.grid = element_blank()
    , plot.subtitle = ggtext::element_markdown(size = 10)
    , strip.text = element_text(color = "black", face = "bold", size = 10)
  )

2.7.3 Objective vs. Treatable Plot

obj_trtbl_temp <-
  constrained_by_scnro_ls_cnstrnt_temp %>% 
  dplyr::filter(cnstrnt_class!="high constraint") %>% 
  dplyr::left_join(
    # aggregate by area_name, scenario_lab
    constrained_by_scnro_ls_cnstrnt_temp %>% 
        dplyr::filter(cnstrnt_class!="high constraint") %>% 
        dplyr::group_by(area_name, scenario_id) %>% 
        dplyr::mutate(low=ifelse(cnstrnt_class=="low constraint",pct_treatable,0)) %>% 
        dplyr::summarise(
          pct_treatable_lowmed = sum(pct_treatable)
          , pct_treatable_low = sum(low)
        ) %>% 
        dplyr::mutate(
          treatment_objective = dplyr::case_when(
              pct_treatable_low >= .2 ~ 1
              , pct_treatable_lowmed >= .2 ~ 2
              , TRUE ~ 3
            ) %>% 
            factor(
              levels = 1:3
              , labels = c("mechanical alone","mechanical+fire","fire")
              , ordered = T
            )
        ) %>% 
        dplyr::group_by(scenario_id) %>% 
        dplyr::arrange(scenario_id, treatment_objective
                       , desc(pct_treatable_low), desc(pct_treatable_lowmed)
        ) %>% 
        dplyr::mutate(sorting_rank = dplyr::row_number()) %>% 
        dplyr::ungroup()
      , by = dplyr::join_by(area_name, scenario_id)
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(
    dplyr::across(
      c(cnstrnt_class, scenario_lab)
      , ~ forcats::fct_rev(.x)
    )
  ) 
# plot
obj_plt_fn <- function(scnum){
plt <- ggplot(
    data = obj_trtbl_temp %>% dplyr::filter(scenario_id==scnum)
    , mapping = aes(
        x = pct_treatable
        , y = reorder(area_name, desc(sorting_rank))
      )
  ) +
  geom_rect(mapping = aes(xmin = 0.2, xmax=0.4, ymin = -Inf, ymax = Inf)
    , fill = "gray77"
    , alpha = 0.6
  ) +
  annotate("text"
     , x = 0.46
     , y = obj_trtbl_temp %>% 
          dplyr::filter(scenario_id == scnum & sorting_rank==max(obj_trtbl_temp$sorting_rank)) %>% 
          dplyr::pull(area_name) %>% 
          unique()
     , label = expression(bold("objective")~bold("20-40%"))
     , size = 3, color = "gray65", parse = T
    ) +
  geom_col(
    mapping = aes(fill = cnstrnt_class)
    , color = NA, width = 0.7
  ) +
  geom_text(
    mapping = aes(
      group = cnstrnt_class
      , label = scales::percent(
          ifelse(pct_treatable<.06,NA,pct_treatable)
          , accuracy = 1
        )
      # , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.2
  ) +
  # total
  geom_text(
    data = obj_trtbl_temp %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::group_by(area_name, scenario_id) %>% 
        dplyr::filter(dplyr::row_number()==1) %>% 
        dplyr::ungroup()
    , mapping = aes(
      x = pct_treatable_lowmed
      , y = reorder(area_name, desc(sorting_rank))
      , label = scales::percent(pct_treatable_lowmed, accuracy = 1) 
      , color = treatment_objective
      , fontface = "bold"
    )
    , size = 3.5
    , hjust = -0.1
    , show.legend = F
  ) +
  # facet_grid(
  #   rows = vars(scenario_lab)
  #   , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
  #   , scales = "free_y"
  #   # , switch = "y"
  # ) +
  scale_fill_manual(values = RColorBrewer::brewer.pal(n=3,name="RdYlBu")[2:3]) +
  scale_color_manual(values = c(
    "mechanical alone"="darkcyan"
    ,"mechanical+fire"="gray33"
    , "fire"="orangered3")
  ) +
  # scale_x_continuous(expand = expansion(mult = c(0, .15)),labels = scales::percent_format()) +
  scale_x_continuous(limits = c(0, .78),labels = scales::percent_format(),position = "top") +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "Percent of Total High-Risk Area"
    , subtitle = obj_trtbl_temp %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::pull(scenario_lab) %>% unique()
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=8)
    , axis.text.x = element_blank()
    , axis.ticks.x = element_blank()
    , axis.text = element_text(size = 8)
    , panel.grid.major.x = element_blank()
    , panel.grid.minor.x = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 10)
    , strip.background = element_rect(fill = "gray88")
    # , strip.placement = "outside"
    , plot.margin = margin(0, 0, 0, 0, "cm")
    , plot.subtitle = element_text(vjust = -4,face="bold")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )
return(plt)
}
cowplot::plot_grid(
  obj_plt_fn(1) +
    theme(
      legend.position = c(0.86,0.82)
      , legend.direction = "vertical"
    )
  , obj_plt_fn(2) +
    theme(legend.position = "none")
  , obj_plt_fn(3) +
    # scale_x_continuous(position = "top") +
    # labs(x = "Percent of Total High-Risk Area") +
    theme(legend.position = "none", plot.margin = margin(0, 0, 0.1, 0, "cm"))
  , ncol = 1
  , align = 'vh'
)

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_obj_sc.png")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

2.8 Totals

##################
# Data Prep
##################
# union landscape and pa level data
# aggregate by scenario for hi-risk pa's
pct_rdctn_pas_temp <- 
  constrained_by_scnro_ls_pa %>%
  dplyr::group_by(pa_id,scenario_id) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::select(scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
                , (tidyselect::starts_with("rmn") & tidyselect::ends_with("area_ha"))
                , covertype_area_ha, feature_area_ha
  ) %>%
  dplyr::group_by(scenario_id, scenario_desc, scenario_lab) %>% 
  dplyr::summarise(
    dplyr::across(
      tidyselect::ends_with("area_ha")
      , sum
    )
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("rmn")
    , names_to = "constraint"
    , values_to = "rmn_area_ha"
    , names_prefix = "rmn"
    , values_drop_na = F
  ) %>%
  dplyr::mutate(constraint = stringr::str_remove(constraint,"_area_ha")) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , pct_rmn = rmn_area_ha/covertype_area_ha
    , pct_rdctn = dplyr::case_when(
      constraint_lvl == 1 ~ -1*(1-pct_rmn)
      , TRUE ~ -1*(dplyr::lag(pct_rmn)-pct_rmn)
    )
    , pct_rdctn_total = max(-1*(1 - ifelse(constraint_lvl == 5, pct_rmn,0)))
  ) %>% 
  dplyr::ungroup()
# overall
pct_rdctn_ls_temp <- 
  constrained_by_scnro_ls %>%
  dplyr::select(scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
                , (tidyselect::starts_with("rmn") & tidyselect::ends_with("area_ha"))
                , covertype_area_ha, feature_area_ha
  ) %>%
  dplyr::group_by(scenario_id, scenario_desc, scenario_lab) %>% 
  dplyr::summarise(
    dplyr::across(
      tidyselect::ends_with("area_ha")
      , sum
    )
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("rmn")
    , names_to = "constraint"
    , values_to = "rmn_area_ha"
    , names_prefix = "rmn"
    , values_drop_na = F
  ) %>%
  dplyr::mutate(constraint = stringr::str_remove(constraint,"_area_ha")) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , pct_rmn = rmn_area_ha/covertype_area_ha
    , pct_rdctn = dplyr::case_when(
      constraint_lvl == 1 ~ -1*(1-pct_rmn)
      , TRUE ~ -1*(dplyr::lag(pct_rmn)-pct_rmn)
    )
    , pct_rdctn_total = max(-1*(1 - ifelse(constraint_lvl == 5, pct_rmn,0)))
  ) %>% 
  dplyr::ungroup()
# combine
pct_rdctn_all_temp <- 
  pct_rdctn_pas_temp %>% 
  dplyr::mutate(dta_src = "High-Risk Project Areas") %>% 
    dplyr::bind_rows(
      pct_rdctn_ls_temp %>% 
        dplyr::mutate(dta_src = "Overall Landscape Area")
    ) %>% 
    dplyr::mutate(
      constraint_lab = factor(
          constraint_lvl
          , ordered = TRUE
          , 1:5
          , labels = c(
            "Protected"
            , "Slope\nSteepness"
            , "Road\nDistance"
            , "Riparian\nBuffer"
            , "Administrative"
          )
        ) %>% forcats::fct_rev()
    ) %>% 
    dplyr::relocate(dta_src)

2.8.1 Reduction Treatable Area Table

tot_table_temp <- pct_rdctn_all_temp %>% 
  dplyr::mutate(constraint_nm = paste0(constraint_lvl,constraint)) %>% 
  dplyr::select(-c(constraint_lvl,constraint,constraint_lab)) %>% 
  tidyr::pivot_wider(
    names_from = constraint_nm
    , values_from = c(rmn_area_ha, pct_rmn, pct_rdctn)
  ) %>% 
  dplyr::mutate(
    pct_covertype_area = covertype_area_ha/feature_area_ha
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::contains("area_ha")
      , ~ scales::comma(.x,suffix = " M", scale = 1e-6, accuracy = 0.1)
    )
    , dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ scales::percent(.x, accuracy = 0.1)
    )
  ) %>% 
  dplyr::arrange(
    desc(dta_src)
    , scenario_id
  )
# make table 
  kableExtra::kable(
      tot_table_temp %>% 
        dplyr::select(
          scenario_lab
          , covertype_area_ha
          , pct_covertype_area
          , pct_rdctn_1protected
          , pct_rdctn_2slope
          , pct_rdctn_3roads
          , pct_rdctn_4riparian
          , pct_rdctn_5administrative
          , rmn_area_ha_5administrative
          , pct_rmn_5administrative
        )
      , caption = paste0("<b>Wildfire Crisis Strategy Priority Landscapes and High-Risk Fireshed Project Areas</b><br>percent reduction of different types of constraints on mechanical treatment")
      , col.names = c(
        ""
        , "Forest+Shrub\n(ha)", "Forest+Shrub\n%"
        , "Protected"
        , "Slope\nSteepness"
        , "Road\nDistance"
        , "Riparian\nBuffer"
        , "Administrative"
        , "Remaining (ha)"
        , "Remaining (%)"
      )
      , escape = F
    ) %>%  
    add_header_above(c(" " = 1, "Strategy Priority\nArea"=2, "Constraints on\nMechanical Treatment" = 5, " " = 2)) %>%
    kable_classic(full_width=T) %>% 
    pack_rows(index = table(forcats::fct_inorder(tot_table_temp$dta_src))) %>% 
    kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%  
    kableExtra::scroll_box(width = "740px")

Table 2.4: Wildfire Crisis Strategy Priority Landscapes and High-Risk Fireshed Project Areas
percent reduction of different types of constraints on mechanical treatment

	Strategy Priority Area		Constraints on Mechanical Treatment
	Forest+Shrub (ha)	Forest+Shrub %	Protected	Slope Steepness	Road Distance	Riparian Buffer	Administrative	Remaining (ha)	Remaining (%)
Overall Landscape Area
Scenario 1	15.8 M	83.1%	-21.9%	-16.9%	-19.9%	-7.0%	-5.5%	4.6 M	28.9%
Scenario 2	15.8 M	83.1%	-21.9%	-5.2%	-12.1%	-10.5%	-8.2%	6.7 M	42.2%
Scenario 3	15.8 M	83.1%	-21.9%	-5.2%	-12.1%	-7.5%	0.0%	8.4 M	53.4%
High-Risk Project Areas
Scenario 1	8.7 M	82.5%	-17.5%	-21.1%	-18.2%	-8.5%	-6.3%	2.5 M	28.5%
Scenario 2	8.7 M	82.5%	-17.5%	-6.4%	-10.5%	-12.9%	-9.4%	3.7 M	43.2%
Scenario 3	8.7 M	82.5%	-17.5%	-6.4%	-10.5%	-9.3%	0.0%	4.9 M	56.3%

2.8.2 Reduction by Constraint

# plot
# plot
ggplot(data = pct_rdctn_all_temp) +
  geom_col(
    mapping = aes(y = scenario_lab, x = pct_rdctn, fill = constraint_lab)
    , color = NA, width = 0.6
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = scenario_lab, x = pct_rdctn
      , label = scales::percent(ifelse(pct_rdctn<.06*-1,pct_rdctn,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = pct_rdctn_all_temp %>% dplyr::filter(constraint_lvl==5)
    , mapping = aes(
      y = scenario_lab, x = pct_rdctn_total
      , label = scales::percent(pct_rdctn_total, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_grid(
      cols = vars(dta_src %>% forcats::fct_rev())
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    ) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0, .12)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_blank()
    , strip.text.x = element_text(color = "black", face = "bold")
    , strip.placement = "outside"
    
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

2.8.3 Distribution of Fireshed PA Constraint

cnstrnt_dta_temp <- constrained_by_scnro_ls_pa %>% 
  dplyr::group_by(pa_id,scenario_id) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::filter(fireshed_crisis_strategy  %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::count(scenario_id, scenario_desc, scenario_lab, cnstrnt_class) %>% 
  dplyr::group_by(scenario_desc, scenario_lab) %>% 
  dplyr::mutate(
    pct = n/sum(n)
    , dta_src = "High-Risk Firesheds"
  ) %>% 
  dplyr::bind_rows(
    constrained_by_scnro_ls_pa %>% 
      dplyr::group_by(pa_id,scenario_id) %>% 
      dplyr::filter(dplyr::row_number()==1) %>% 
      dplyr::ungroup() %>% 
      dplyr::count(scenario_id, scenario_desc, scenario_lab, cnstrnt_class) %>% 
      dplyr::group_by(scenario_desc, scenario_lab) %>% 
      dplyr::mutate(
        pct = n/sum(n)
        , dta_src = "Overall Landscape Area"
      )
  ) %>% 
  dplyr::mutate(
    dta_src = dta_src %>% forcats::fct_rev()
  )
# n for notation
nlab_temp <- 
  cnstrnt_dta_temp %>% 
  dplyr::filter(scenario_id==1) %>% 
  dplyr::group_by(dta_src) %>% 
  dplyr::summarise(n = sum(n)) %>% 
  dplyr::mutate(l = paste0(dta_src, " n = ", scales::comma(n))) %>% 
  dplyr::pull(l) %>% 
  paste(collapse = "\n")
# plot base
plt_fshed_cnstrnt_lvl_temp <- ggplot(data=cnstrnt_dta_temp) +
  geom_col(
    mapping = aes(x = pct, y = scenario_lab, fill=cnstrnt_class)
    ,width = 0.7, alpha=0.8
  ) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  facet_grid(cols = vars(dta_src)) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "% of Fireshed\nProject Areas"
    # , caption = nlab_temp
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
    , plot.caption = element_text(size = 8)
  ) + 
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )
# plot with large labels
plt_fshed_cnstrnt_lvl_temp + 
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.065,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 4 # 2.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = ifelse(pct>=0.065,
          paste0("n="
            ,scales::comma(n, accuracy = 1)
          )
          , NA)
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.5 # 2.8
    , vjust = 2.3
  ) 

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_prjareas_dist_sc.png")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

2.8.3.1 Combine map and fireshed plot

# plot small for inset
inset_temp <- plt_fshed_cnstrnt_lvl_temp +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.105,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3.2 # 2.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = ifelse(pct>=0.105,
          paste0("n="
            ,scales::comma(n, accuracy = 1)
          )
          , NA)
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.1 # 2.8
    , vjust = 2.3
  ) +
    cowplot::theme_minimal_grid(9) +
    theme(
      legend.position = "none"
      , plot.background = element_rect(color = "black", fill="white", linewidth = 2)
      , axis.text.x = element_blank()
      , axis.title.x = element_text(size=7)
      , strip.text = element_text(color = "black", face = "bold", size = 9)
    )

# combine
cowplot::ggdraw(
  plt_fpas_sc_map[[2]] + 
    theme(
      legend.position = c(.4,.98)
      , legend.direction = "horizontal"
      , legend.text = element_text(size = 9)
      , plot.background = element_rect(fill="white", color = NA)
      , plot.margin = margin(0, 5.1, -0.7, 0.1, "cm")
    )
  ) +
  cowplot::draw_plot(inset_temp, x=.5, y=.56, width=.45, height=.4) +
  cowplot::draw_plot_label(
    label = c("A", "B"),
    x = c(0.02, 0.5),
    y = c(1, 0.95),
    size = 17,
    fontface = "bold",
    color = "gray11"
  )

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_prjareas_sc_map_inset.png")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

2.9 Landcover Class Area

2.9.1 Overall Landscape

What is the landcover classification of the area within the landscape boundary? This GEE script was used to calculate the area by NLCD landcover class.

# nlcd class groups
nlcd_reclass_df_temp <- data.frame(
    class_id = c(11,12,21,22,23,24,31,41,42,43,51,52,71,72,73,74,81,82,90,95) %>% as.character()
    , class_label = c(
      "Water"
      , "Water"
      , "Developed"
      , "Developed"
      , "Developed"
      , "Developed"
      , "Barren"
      , "Forest"
      , "Forest"
      , "Forest"
      , "Shrubland"
      , "Shrubland"
      , "Herbaceous"
      , "Herbaceous"
      , "Herbaceous"
      , "Herbaceous"
      , "Planted/Cultivated" # double-check name
      , "Planted/Cultivated"
      , "Wetlands"
      , "Wetlands"
      )
  ) %>% 
  dplyr::left_join(
    FedData::nlcd_colors() %>% dplyr::rename(class_id=ID) %>% dplyr::rename_with(tolower) %>% dplyr::mutate_all(as.character)
    , by = dplyr::join_by(class_id)
  )
# load data by wf priority 
wf_nlcd_area <- readr::read_csv("../data/wildfirepriority/wfpriority_nlcd_area.csv") %>% 
  dplyr::rename_with(tolower) %>% 
  dplyr::left_join(
    wf_landscapes %>% 
      dplyr::select(name,state,area_name)
    , by = dplyr::join_by(name,state)
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("area_m2_nlcd_cl")
    , names_to = "class_id"
    , values_to = "area_m2"
    , names_prefix = "area_m2_nlcd_cl_"
    , values_drop_na = F
  ) %>% 
  dplyr::left_join(
    nlcd_reclass_df_temp
    , by = dplyr::join_by(class_id)
  ) %>% 
  dplyr::mutate(
    class_label = forcats::fct_relevel(class_label, c("Forest","Shrubland")) %>% forcats::fct_rev()
  )
# get color palette
nlcd_colors <- wf_nlcd_area %>% 
  dplyr::group_by(class_id,color,class_label) %>% 
  dplyr::summarise(sum_area = sum(area_m2,na.rm = T)) %>% 
  dplyr::arrange(class_label,desc(sum_area)) %>% 
  dplyr::group_by(class_label) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::pull(color)
# aggregate and plot
wf_nlcd_area %>% 
  dplyr::group_by(area_name,class_label) %>% 
  dplyr::summarise(sum_area_ha = sum(area_m2,na.rm = F)/10000) %>% 
  dplyr::mutate(sum_area_ha = ifelse(is.na(sum_area_ha),0,sum_area_ha)) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(
    pct_area = sum_area_ha/sum(sum_area_ha)
    # , area_name = area_name %>% forcats::fct_rev()
  ) %>%
  dplyr::ungroup() %>% 
ggplot(
    mapping = aes(x = pct_area, y = reorder(area_name, desc(area_name)), fill=class_label)
  ) +
  geom_col(width = 0.7, alpha=0.8) +
  geom_text(
    mapping = aes(
      label = scales::percent(ifelse(pct_area>.03,pct_area,NA), accuracy = 1)
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  scale_fill_manual(values = nlcd_colors) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "NLCD\nCover\nType"
    , y = ""
    , x = "% Landscape Area"
  ) +
  theme_light() +
  theme(
    legend.position = "top" # c(0.9, 0.9)
    , legend.title = element_text(size=9)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

2.9.2 Fireshed Project Areas

What is the landcover classification of the fireshed project area units? This GEE script was used to calculate the area by NLCD landcover class.

# load nlcd cover data by fireshed project area
nlcd_area_ls_pa <- 
  readr::read_csv("../data/wildfirepriority/wfpriority_fireshed_nlcd_area.csv") %>% 
  dplyr::rename_with(tolower) %>% 
  dplyr::select(
    pa_id
    , tidyselect::contains("area_m2")
  ) %>% 
  dplyr::mutate(pa_id = as.character(pa_id)) %>% 
  dplyr::inner_join(
    constrained_by_scnro_ls_pa %>% 
      dplyr::filter(scenario_id==1) %>% 
      dplyr::select(pa_id,name,state,area_name,fireshed_crisis_strategy)
    , by = dplyr::join_by(pa_id)
    , multiple = "all"
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("area_m2_nlcd_cl")
    , names_to = "class_id"
    , values_to = "area_m2"
    , names_prefix = "area_m2_nlcd_cl_"
    , values_drop_na = F
  ) %>% 
  dplyr::left_join(
    nlcd_reclass_df_temp
    , by = dplyr::join_by(class_id)
  ) %>% 
  dplyr::mutate(
    class_label = forcats::fct_relevel(class_label, c("Forest","Shrubland")) %>% forcats::fct_rev()
  )

# aggregate and plot
nlcd_area_ls_pa %>% 
  dplyr::group_by(area_name,class_label) %>% 
  dplyr::summarise(sum_area_ha = sum(area_m2,na.rm = F)/10000) %>% 
  dplyr::mutate(sum_area_ha = ifelse(is.na(sum_area_ha),0,sum_area_ha)) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(pct_area = sum_area_ha/sum(sum_area_ha)) %>%
  dplyr::ungroup() %>% 
ggplot(
    mapping = aes(x = pct_area, y = reorder(area_name, desc(area_name)), fill=class_label)
  ) +
  geom_col(width = 0.7, alpha=0.8) +
  geom_text(
    mapping = aes(
      label = scales::percent(ifelse(pct_area>.03,pct_area,NA), accuracy = 1)
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  scale_fill_manual(values = nlcd_colors) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "NLCD\nCover\nType"
    , y = ""
    , x = "% Landscape Area"
    , subtitle = "All Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top" # c(0.9, 0.9)
    , legend.title = element_text(size=9)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

The figure above shows the NLCD landcover class distribution for all fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape.

2.9.3 NLCD Cover Type high-risk firesheds

# aggregate and plot
nlcd_area_ls_pa %>%
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::group_by(area_name,class_label) %>% 
  dplyr::summarise(sum_area_ha = sum(area_m2,na.rm = F)/10000) %>% 
  dplyr::mutate(sum_area_ha = ifelse(is.na(sum_area_ha),0,sum_area_ha)) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(pct_area = sum_area_ha/sum(sum_area_ha)) %>%
  dplyr::ungroup() %>% 
ggplot(
    mapping = aes(x = pct_area, y = reorder(area_name, desc(area_name)), fill=class_label)
  ) +
  geom_col(width = 0.7, alpha=0.8) +
  geom_text(
    mapping = aes(
      label = scales::percent(ifelse(pct_area>.03,pct_area,NA), accuracy = 1)
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  scale_fill_manual(values = nlcd_colors) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "NLCD\nCover\nType"
    , y = ""
    , x = "% Landscape Area"
    , subtitle = "High-Risk Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top" # c(0.9, 0.9)
    , legend.title = element_text(size=9)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

The figure above shows the NLCD landcover class distribution for only high-risk fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape.

2.9.4 NLCD Cover Type vs. Mechanically Available

Is there a relationship between forest and shrubland cover and mechanical availability and is that relationship the same?

nlcd_area_ls_pa_cl_temp <-
  nlcd_area_ls_pa %>% 
    # dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
    dplyr::group_by(pa_id,area_name,class_label) %>% 
    dplyr::summarise(sum_area_ha = sum(area_m2,na.rm = F)/10000) %>% 
    dplyr::group_by(pa_id,area_name) %>% 
    dplyr::mutate(pct_area = sum_area_ha/sum(sum_area_ha)) %>%
    dplyr::ungroup() %>% 
    dplyr::filter(class_label %in% c("Forest","Shrubland")) %>% 
    dplyr::inner_join(
      constrained_by_scnro_ls_pa %>% 
        dplyr::select(
          scenario_id
          , scenario_lab
          , scenario_desc
          , pa_id
          , area_name
          , pct_rmn5_administrative
          , pct_rdctn_total
          , cnstrnt_class
          , fireshed_crisis_strategy
        )
      , by = dplyr::join_by(pa_id,area_name)
      , multiple = "all"
    )  
  # dplyr::mutate(scenario_desc = scenario_desc %>% forcats::fct_rev())
# treatable vs covertype
ggplot(
    data = nlcd_area_ls_pa_cl_temp
    , mapping = aes(y = pct_rmn5_administrative, x = pct_area, color = class_label)
  ) +
  geom_point(alpha = 0.7, size = 1) +
  # geom_smooth(method = "lm", lwd = 1, color = "gray", linetype = "dashed", se=F) +
  scale_color_manual(values = nlcd_colors[(length(nlcd_colors)-1):length(nlcd_colors)]) +
  facet_grid(cols = vars(class_label), rows = vars(scenario_desc)) +
  scale_x_continuous(labels = scales::percent_format()) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    x = "Percent Cover Type"
    , y = "Percent Treatable Area Remaining"
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , axis.title = element_text(size=9)
    , axis.text = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
  )

All fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape shown in figure above.

2.9.5 Distribution of Mechanically Available (Shrb vs. Frst)

# n for notation
nlab_temp <- nlcd_area_ls_pa_cl_temp %>% 
  dplyr::filter(pct_area > .5 & scenario_id==1) %>% 
  dplyr::count(class_label) %>% 
  dplyr::mutate(l = paste0(class_label, " n = ", scales::comma(n))) %>% 
  dplyr::pull(l) %>% 
  paste(collapse = " | ")
# plot
nlcd_area_ls_pa_cl_temp %>% 
  dplyr::filter(pct_area > .5) %>% 
ggplot(
    mapping = aes(x = pct_rmn5_administrative, fill = class_label, group = class_label, color = class_label)
  ) +
  geom_density(alpha = 0.7) +
  scale_fill_manual(values = nlcd_colors[(length(nlcd_colors)-1):length(nlcd_colors)]) +
  scale_color_manual(values = nlcd_colors[(length(nlcd_colors)-1):length(nlcd_colors)]) +
  facet_grid(cols = vars(class_label), rows = vars(scenario_desc)) +
  scale_x_continuous(labels = scales::percent_format()) +
  # scale_y_continuous(labels = scales::percent_format()) +
  labs(
    y = "Density"
    , x = "Percent Treatable Area Remaining"
    , caption = nlab_temp
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , axis.title = element_text(size=9)
    , axis.text = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
    , plot.caption = element_text(size = 8)
  )

In the plot above fireshed project areas are grouped into forest or shrubland based on the majority cover type (e.g. >50% forest). Fireshed project areas where the majority covertype is not forest or shrubland are not shown. All fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape shown in figure above.

2.9.6 Distribution of Fireshed PA Constraint (Shrb vs. Frst)

nlcd_area_ls_pa_cl_temp %>% 
  dplyr::filter(pct_area > .5) %>% 
  dplyr::count(class_label, scenario_desc, scenario_lab, cnstrnt_class) %>% 
  dplyr::group_by(class_label, scenario_desc, scenario_lab) %>% 
  dplyr::mutate(pct = n/sum(n)) %>% 
ggplot() +
  geom_col(
    mapping = aes(x = pct, y = scenario_lab, fill=cnstrnt_class)
    ,width = 0.7, alpha=0.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.05,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  facet_grid(cols = vars(class_label)) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "% of Fireshed\nProject Areas"
    , caption = nlab_temp
    , subtitle = "All Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
    , plot.caption = element_text(size = 8)
  ) + 
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

In the plot above fireshed project areas are grouped into forest or shrubland based on the majority cover type (e.g. >50% forest). Fireshed project areas where the majority covertype is not forest or shrubland are not shown. All fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape shown in figure above.

2.9.7 High-Risk Distribution of Fireshed PA Constraint (Shrb vs. Frst)

# n for notation
nlab_temp <- nlcd_area_ls_pa_cl_temp %>% 
  dplyr::filter(pct_area > .5 
                & fireshed_crisis_strategy  %in% c("USFS-Only", "All-Lands") 
                & scenario_id==1) %>% 
  dplyr::count(class_label) %>% 
  dplyr::mutate(l = paste0(class_label, " n = ", scales::comma(n))) %>% 
  dplyr::pull(l) %>% 
  paste(collapse = " | ")
# plot
nlcd_area_ls_pa_cl_temp %>% 
  dplyr::filter(pct_area > .5 & fireshed_crisis_strategy  %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::count(class_label, scenario_desc, scenario_lab, cnstrnt_class) %>% 
  dplyr::group_by(class_label, scenario_desc, scenario_lab) %>% 
  dplyr::mutate(pct = n/sum(n)) %>% 
ggplot() +
  geom_col(
    mapping = aes(x = pct, y = scenario_lab, fill=cnstrnt_class)
    ,width = 0.7, alpha=0.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.05,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  facet_grid(cols = vars(class_label)) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "% of Fireshed\nProject Areas"
    , caption = nlab_temp
    , subtitle = "High-Risk Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
    , plot.caption = element_text(size = 8)
  ) + 
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

In the plot above fireshed project areas are grouped into forest or shrubland based on the majority cover type (e.g. >50% forest). Fireshed project areas where the majority covertype is not forest or shrubland are not shown. Only high-risk fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape shown in figure above.